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ABSTRACT 


Cave  systems  with  two  cr  more  entrances  at  different 
elevations  often  exhibit  seasonally  directional  air  flow. 
Tnis  flow  is  analogous  to  that  in  a  chimney  and  is  thus 
referred  to  as  the  "chimney”  or  "stack"  effect.  Previous 
attempts  tc  model  cave  winds  of  this  type  have  greatly 
simplified  the  phenomenon  by  assuming  static  conditions  and 
then  estimating  air  pressure  or  density  differences  between 
the  cave  atmosphere  and  the  outside  air.  This  approach  is 
basically  incorrect  when  applied  to  dynamic  conditions,  and 
may  lead  to  substantial  errors  in  estimating 
speleomicroclimatic  parameters.  Wigley  and  Brown  (1972)  have 
predicted  temperature  profiles  under  dynamic  conditions,  and 
this  thesis  extends  their  work  to  a  general  model  for  this 
type  of  air  flow.  The  assumption  of  the  incompressibility  of 
air  is  analysed  and  found  to  be  a  reasonable  simplification. 
The  results  are  of  significance  to  geomorphologists 
investigating  inaccessible  karst  systems,  biospeleologists 
concerned  with  the  distribution  of  life  in  caves,  and  mine 
ventilation  engineers. 
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I  -  CAVES  AND  THE  CAVE  ATMOSPHEEE 


From  very  early  in  human  prehistory,  man  has  been 
associated  with  caves.  They  have  provided  him  with  shelter 
from  the  elements,  protection  from  his  enemies,  and  often  a 
means  of  living  by  the  extraction  of  valuable  minerals.  Even 
in  the  heavily  industrialised  western  world,  they  are  still 
used  as  storage  areas  for  perishables  such  as  wine,  and  as  a 
source  of  tat  guano,  an  excellent  organic  fertilizer.  More 
recently,  however,  caves  have  ccme  under  the  scrutiny  of  the 
scientific  community.  Aside  from  the  obvious  concerns  of  how 
caves  are  formed  and  what  effect  they  have  on  local  and 
regional  drainage  patterns,  an  aspect  of  speleology  that  has 
elicited  much  interest  is  spelecmicroclimatology,  the  study 
of  the  cave  atmosphere.  Of  considerable  interest  is  the 
relative  invariability  of  cave  atmospheres,  since  this 
invariability  is  one  of  the  most  useful  features  of  caves. 
Because  internal  temperatures  approximate  the  mean  annual 
temperature  outside  a  cave,1  caves  which  are  located  in 
regions  of  significant  seasonal  temperature  variations  are 
substantiall j  warmer  than  the  exterior  region  in  winter  and 
substantially  cooler  in  summer.  As  a  means  of  surviving  the 
harsh  winters  of  early  post-glacial  times,  for  example, 
caves  were  nearly  ideal. 


i  As  will  be  seen  later,  internal  temperatures  may  vary 
significantly  near  entrances  which  draw  air,  but  the 
amplitude  of  these  fluctuations  decreases  away  from  an 
entrance  to  the  point  where  temperatures  are  nearly 
invariant. 
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The  cave  atmosphere  also  determines,  to  a  large  extent, 
the  flora  and  fauna  found  beneath  the  surface,  Poulson  and 
White  (1969)  have  characterised  caves  as  ecological  and 
evolutionary  laboratories  where  relations  between  surface 
and  subsurface  species,  as  well  as  relations  within  sub¬ 
surface  communities,  may  be  more  easily  studied  than  those 
amongst  surface  species  and  communities  (see  also  Vandel, 
1965)  . 

Karst  hydrologists  have  for  many  years  used  the  flow  of 
groundwater  to  study  inaccessible  karst  aquifers.  The 
rapidity  with  which  water  and  water  waves  flow  through  the 
aquifer,  as  well  as  the  temperature  and  geochemical 
characteristics  of  the  water  entering  and  leaving  the 
aquifer,  yield  information  concerning  the  character  of  that 
aquifer  (see  Erown  et  al.,  1969;  Brown  1972;  Shuster  and 
White,  1971;  among  others).  The  possibility  of  gleaning 
similar  information  from  the  flow  of  air  has  been  little 
studied,  due  in  part  to  the  substantially  increased 
complexity  of  the  mathematical  description  of  the  flow. 
Wigley  (Wigley  et  al. ,  1966;  Wigley,  1967;  Wigley  and  Brown, 
1969)  has  done  considerable  work  on  estimating  total 
permeability  of  aquifers  from  velocity  and  pressure 
variations  of  air  at  the  entrance  to  a  cave  in  Australia, 
while  Plummer  (1962,  1969)  has  attempted  estimates  of  cave 
volume  from  the  resonance  of  certain  caves.  Both  these 
studies,  however,  deal  with  non-dynamic  caves,  and  work  on 
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dynamic  caves  is  virtually  non-existent.  Since  a  knowledge 
of  the  microclimatic  conditions  to  be  expected  in  dynamic 
caves  would  be  vital  to  the  use  of  air-flow  in  the  remote 
sensing  of  karst  aquifer  characteristics,  a  dynamic  model 
should  be  of  interest  to  karst  hydrologists. 

Much  akin  to  the  cave  atmosphere  is  the  mine 
atmosphere,  which  is  of  great  interest  to  the  mining 
engineer.  To  ensure  safe  and  comfortable  working  conditions, 
the  natural  microclimatic  conditions  found  in  mines  must  be 
ascertained  and,  if  necessary,  altered  by  mechanical 
ventilation.  While  mines  are  found  in  a  much  broader 
spectrum  of  climatic  and  geologic  conditions,  the  same  basic 
principles  are  responsible  for  the  determination  of 
microclimatic  conditions  in  both  mines  and  caves. 

Basic  Principles 

Having  indicated  the  importance  of  subterranean 
microclimatic  conditions,  one  might  expect  a  question  as  to 
what  determines  these  conditions.  The  answer  has  been 
intimated  by  the  statement  that  cave  temperatures 
approximate  the  mean  annual  surface  temperature;  hence  the 
cave  atmosphere  is  not  a  closed  system,  but  rather  is 
directly  connected  with  and  to  a  large  extent  determined  by 
the  free  atmosphere.  The  nature  of  this  relationship  is 
dependent  upon  the  morphology  of  the  given  cave,  especially 
the  number  of  entrances-  The  classical  division  of  caves 


. 


■  ■ 

- 


. 

:  IS  J 


4 


into  "static"  and  "dynamic"  was  popularised  by  Geiger  (see 
for  instance  Geiger ,  1966).  By  implication,  static  caves 
(those  with  only  one  entrance)  do  not  vary  microclimatic¬ 
ally.  Wigley  and  Brown  (in  press)  point  out  that  this 
definition  is  inadequate,  for  one-entrance  cave  systems 
often  have  very  active  air  circulations  (see  Conn,  1966; 
Wigley  et  al.,  1966;  Wigley,  1967;  Wigley  and  Brown,  1969) , 
but  for  simplicity,  this  terminology  will  be  maintained. 
Figure  1  shows  an  idealised  static  cave.  As  external 
pressure  (P[ext])  rises,  the  internal  pressure  measured  at 
the  surface  (P£int])  must  also  rise  to  maintain  equilibrium. 
Otherwise  a  pressure  differential  could  exist  at  the 
entrance.  In  crder  for  P£int]  to  increase,  air  must  flow 
into  the  cave  from  the  surface.  This  air  circulation  is 
fundamental  to  the  determination  of  speleomicroclimatic 
conditions:  "Wind  is  not  only  the  most  striking  phenomenon 
of  weather  in  caves,  it  is  also  its  determining  factor..." 
(Bogli  and  Franke,  1968,  p.  39).  Likewise,  if  P£ext  ]  < 
P£int],  there  is  a  flow  of  air  out  the  entrance.  Thus  air 
circulation  and  its  effect  on  speleomicroclimatic  conditions 
are,  in  static  caves,  related  primarily  to  external 
barometric  pressure  changes. 

"Dynamic"  caves  are  those  with  two  or  more  entrances 
vertically  separated  by  at  least  a  few  meters.  A  markedly 
different  relationship  between  the  cave  and  free  atmosphere 
exists  here,  for  air  need  not  flow  in  and  out  the  same 
entrance.  As  previously  stated,  speleomicroclimatic 
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parameters  are  considerably  less  variable  than  their  surface 
counterparts.  This  is  due  in  large  part  to  the  relative 
invariability  of  the  rock  temperature  within  the  cave. 
Assume  for  the  moment  that  the  internal  and  external 
atmospheres  are  balanced;  that  is,  they  are  in  static 
equilibrium  (see  Figure  2) .  Now  increase  the  average 
external  temperature.  A  pressure  imbalance  is  created  such 
that  air  flows  into  the  upper  entrance  and  out  the  lower 


entrance.  If  a 

static 

model  is 

assumed. 

the 

pressure 

difference  at 

either 

entrance 

(i  •  e  •  , 

the 

pressure 

differential  that  would  exist  across  any  barrier  within  the 
cave  that  completely  stepped  the  flow  of  air)  is  the 
difference  in  weight  of  a  column  of  cave  air  and  a  column  of 
external  air  whose  heights  are  the  height  difference  between 
the  two  entrances  (Figure  3) .  As  the  warmer  air  flows  into 
the  cave,  it  is  cooled  by  the  rock  walls  of  the  cave  until 
it  is  at  virtually  the  same  temperature  as  the  rock.  At  the 
same  time,  the  rock  walls  are  warmed  by  the  air,  but  this 
process  is  much  slower  because  of  the  different  heat 
capacity  of  the  rock.  If  external  conditions  did  not  again 
vary,  then  the  cave  and  free  atmospheres  would,  over  a 
period  of  years,  once  more  achieve  static  equilibrium. 
External  conditions  obviously  vary  seasonally,  however,  and 
over  such  a  short  period  of  time  there  is  very  little 
variation  in  the  rock  wall  temperature.  Air  circulation, 
therefore,  is  strongly  directional  with  air  flowing  steadily 
into  the  upper  entrance  and  out  the  lower  entrance  when 
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external  temperatures  are  greater  than  the  internal  rock- 
wall  temperature,  a  condition  found  consistently  in  summer. 
In  winter,  external  temperatures  are  lower  than  the  internal 
wall  temperature  and  airflow  is  in  the  opposite  direction. 
Tnis  dynamic  wind  flow  is  called  the  "chimney"  (or  "stack") 
effect  because  of  the  strong  analogy  (under  winter 
conditions)  tc  air  flow  in  a  chimney.  Air  within  the  chimney 
(cave)  is  warmer  than  external  air  so  a  pressure  imbalance 
exists  which  draws  air  into  the  fireplace  (lower  entrance) 
and  pushes  it  out  the  top  of  the  chimney  (upper  entrance) . 

Subject  of  Thesis 

This  thesis  will  be  concerned  specifically  with  the 
theoretical  prediction  of  microclimatic  parameters  in  an 
ideal  dynamic  cave,  given  the  external  atmospheric 
conditions  and  the  cave  morphology.  It  is  recognised  that  in 
practice  the  air  circulation  in  any  particular  cave  may 
often  be  a  combination  of  both  static  and  dynamic 
processes,2  tut  the  concern  here  will  be  for  an  idealised 
two-entrance  system  whose  internal  atmosphere  is  conditioned 
by  dynamic  influences  only.  The  proposed  model  will 
therefore  be  a  dynamic  one  based  on  the  principles  of 
thermodynamics  and  heat  and  mass  transfer.  This  is  a 
significant  departure  from  most  previous  work,  which  has 
assumed  static  or  quasi-static  models  for  the  sake  of 
simplicity.  Beth  compressible  and  incompressible  steady  flow 
will  be  investigated.  Since,  as  previously  stated,  the 
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prediction  of  subterranean  atmospheric  conditions  is  of 
interest  to  speleologists,  biologists,  karst  hydrologists, 
and  mining  engineers  alike,  the  results  of  this  thesis 
should  be  of  considerable  significance  to  a  number  of 
disciplines • 


2  While  barcmetric  pressure  changes  and  the  chimney  effect 
are  the  predominant  causes  of  air  circulation,  other  causes 
should  not  be  overlooked.  Winds  induced  by  the  surface  drag 
of  flowing  water  within  caves  are  very  common  (see  Trombe, 
1952;  Myers,  1962),  and  Helmholtz  resonance  has  been 
observed  in  many  caves  (Plummer,  1962,  1969;  Moore  and 
Nicholas,  1964;  see  also  Schmidt,  1958) .  for  further 
discussion  of  the  types  of  winds  found  in  caves,  see  Trombe 
(1952)  ,  Myers  (1962)  and  Trimmel  (1968)  . 
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PREVIOUS  RESEARCH 


Static  Theories 

Although  the  existence  of  the  chimney  effect  in  caves 
has  been  known  for  years,  few  researchers  have  attempted  to 
describe  it  mathematically.  Most  are  content  to  point  out 
that,  in  summer,  cave  air  is  colder  than  external  air,  so  it 
descends;  while  in  winter,  the  reverse  occurs.  Nearly 
everyone  who  has  attempted  to  model  the  chimney  phenomenon 
has  chosen  to  treat  density  (or  pressure)  differences  in 
terms  of  hydrostatics.  While  it  is  obviously  incorrect  to 
use  hydrostatics  in  a  dynamic  situation,  the  approach  is 
understandable  due  to  a  persistent  lack  of  understanding  of 
the  principles  involved  and  due  to  the  mathematical 
simplicity  which  this  approach  affords. 

At  the  turn  of  the  century,  the  cave  atmosphere  was 
thought  of  as  distinct  from  the  free  atmosphere:  "...static 
influences,  especially  those  due  to  differences  in 
temperature,  were  regarded  as  decisive..."  in  determining 
air  circulation  (Bogli  and  Franke,  1968,  p.  39) .  One  of  the 
first  to  use  this  approach  was  Simonys  (1913) .  He  equated 
the  average  kinetic  energy  of  the  cave  air  with  the  work 
required  to  move  a  parcel  of  air  of  unit  length  through  the 
cave.  This  model,  however,  is  obviously  a  gross 

simplification.  In  his  classical  Traite _ de _ Speleologie, 

Trombe  (1952)  essentially  reiterated  the  conventional  static 
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model.  Bather  than  actually  calculate  velocities,  he  (and 
most  subsequent  proponents  of  the  static  model)  concerned 
himself  with  the  motive  force  behind  the  air  movement, 
deriving  the  formula 

PDIFF  =  (0.1)  (H)  (273)  ( (D£  int  ]/T£  int  ])  -  (D£  ext  ]/T[  ext  ])  ) 

where  PDIFF  is  the  pressure  difference  (in  the  unfortunate 
dimensions  of  gms  force/cm2)  between  internal  and  external 
columns  of  air  (the  motive  force) ;  H  is  the  elevation 
difference  between  the  two  entrances  (meters) , 
D[  int ] (gms/liter)  and  T[int](°K)  are  the  internal  air 
density  (corrected  to  0°C)  and  temperature;  and  D£ext]  and 
T£ext]  are  the  corresponding  external  parameters.  The  D's 
and  T's  are  assumed  to  be  average  values  taken  at  an 
elevation  halfway  between  the  entrances,  on  the  presuppos¬ 
ition  that  they  vary  linearly  with  elevation.  This  is 
definitely  an  over- simplification  in  a  dynamic  situation, 
and  represents  a  basic  flaw  in  Trombe's  approach. 
Nevertheless,  this  very  simplicity  has  attracted  many 
individuals  to  adopt  the  static  model.  Montoriol  Pous  (?)  , 
in  explaining  meteorological  observations  in  a  Spanish  cave, 
adopted  Trombe's  arguments  without  modification.  More 
recently,  however,  Cigna  (1963,  1965)  has  modified  Trombe's 
model  to  include  the  effects  of  humidity  on  air  density. 
Nothing,  however,  is  mentioned  about  density  variations 
other  than  those  directly  attributable  to  changes  in 


elevation. 
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Another  proponent  of  the  static  theory,  Eraso  (1965) , 
has  developed  a  nomogram  for  cave  climate  calculations. 
Basically,  he  treats  cave  and  surface  air  masses  as  two 
separate  entities  which  are  mixed  in  the  entrance  area  to 
produce  an  air  mass  differing  in  temperature  and  density 
from  both  original  air  masses.  The  nomogram  is  intended  to 
simplify  calculations  of  the  resulting  temperature  and 
density.  In  effect,  as  Wigley  and  Brown  (1971)  point  out, 
the  nomogram  states  that  the  transition  or  "cold  zone" 
temperature  (zone  of  significant  temperature  fluctuation 
near  the  surface)  is  equal  to  the  external  adiabatic  wet- 
bulb  temperature,  but  this  is  fundamentally  incorrect: 
"...the  agreement  between  the  theory  of  Eraso  and 
observation  is  [merely]  fortuitous..."  (Wigley  and  Brown, 
1971,  p.  .314)  . 

Dynamic  Theories 

In  the  mid-1960s,  speleologists  began  to  understand 
that  microclimatic  parameters  such  as  temperature  were  not 
constant  within  caves,  but  varied  significantly  even  several 
thousand  feet  from  the  entrance.  The  previous  assumption  was 
that  internal  temperatures  away  from  the  entrance  were 
governed  by  the  transmission  of  surface  conditions  through 
solid  bedrock  (specifically  limestone) .  Since  there  is 
little  effect  from  seasonal  fluctuations  beyond  a  depth  of 
fifty  feet  (see  for  instance  Cropley,  1965) ,3  the  conclusion 
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remained  that  deep  cave  temperatures  are  essentially 
invariable.  While  agreeing  with  this  basic  argument,  Cropley 
(1965)  actually  observed  variations  of  ±1°F  over  a  mile  from 
the  entrance  of  several  caves  in  West  Virginia.  He  concluded 
that  three  temperature  zones  existed  in  these  caves.  Zone  1 
is  the  zone  governed  by  entrance  effects.  Zone  3  is  the  zone 
of  constant  temperature  where  conditions  are  controlled  by 
heat  transmitted  through  limestone  from  the  surface.  Zone  2, 
however,  is  a  dynamic  zone  where  temperatures  vary 
seasonally  with  an  exponentially  decreasing  amplitude  as  one 
progresses  into  the  cave  from  any  entrance  which  draws  air. 

Cropley  does  not  directly  discuss  wind  velocity, 
however,  and  refrains  from  a  theoretical  discussion  of 
temperature  distribution.  But  for  the  purposes  of  this 
thesis,  since  the  chimney  effect  is  dependent  upon  an 
internal-external  temperature  differential,  the  temperature 
distribution  within  a  cave  is  critical  to  the  model  to  be 
developed.  Wigley  and  Brown  (1971)  have  succeeded  in 
adequately  predicting  cave  temperatures  by  including  (in  a 
dynamic  model)  the  effects  of  condensation  and  evaporation, 
based  on  heat  and  mass  transfer  theory.  They  assume  a 
constant  wall  temperature  and  show  that  the  difference 
between  the  air  temperature  and  the  internal  wall 
temperature  decays  quasi-exponentially  to  zero  as  one 

3  Lange  (1954),  extends  this  limit  deeper  than  50  ft., 
the  argument  is  still  the  same. 
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proceeds  along  the  cave  passage  in  the  direction  of  the 
wind.  The  rate  of  decay  was  found  to  be  dependent  upon 
passage  size  and  wind  velocity.  This  study  by  Wigley  and 
Brown  will  serve  as  the  basis  for  the  present  thesis. 


■Wl  $J 
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Ill  -  RELATED  BESEARCH 

A  problem  similar  to  the  determination  of  microclimatic 
parameters  of  the  cave  atmosphere  is  found  in  the  planning 
of  mine  ventilation.  Mines  with  two  or  more  connections  to 
the  surface  are  analogous  to  dynamic  caves.  To  maintain 
healthy  working  conditions,  many  mines  must  be  artificially 
ventilated,  and  the  degree  to  which  the  natural  ventilation 
(due  to  the  chimney  effect)  helps  or  hinders  this  artificial 
ventilation  is  of  considerable  importance. 

Mine  engineers  have  been  concerned  with  mine 
ventilation  for  a  much  longer  period  than  have  micro¬ 
climatologists  with  caves,  and  it  is  only  to  be  expected 
that  their  analysis  of  the  chimney  effect  (and  its  effect 
upon  the  mine  atmosphere)  is  somewhat  more  refined  than  the 
microclimatclogists • .  The  major  advance  of  engineers  in  this 
area  has  been  the  calculation  of  what  is  termed  the  "natural 
ventilation  pressure."  In  the  absence  of  artificially- 
induced  density  differences,  this  may  be  defined  as  the 
pressure  difference  (due  to  density  differences)  between  two 
columns  of  air  of  equal  height,  (technically  the  height  of 
the  atmosphere  above  the  lower  entrance,  but  generally  taken 
as  the  vertical  separation  of  the  two  entrances)  one  within 
and  the  other  exterior  to  the  mine.  This  is  the  shaft-and— 
adit  case,  which  best  approximates  a  dynamic  cave.  In  mine 
ventilation  literature,  however,  the  natural  ventilation 
pressure  (NVP)  also  refers  to  the  pressure  difference 
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resulting  from  density  differences  induced  by  artificial 
ventilation  (as  opposed  to  the  pressure  differential  due  to 
the  fan  itself) .  Because  the  two  entrances  to  the  mine  need 
not  be  at  different  elevations  (see  Figure  4) ,  this 
situation  is  markedly  different  from  dynamic  caves  (where  no 
artificial  ventilation  exists). 

Early  attempts  to  calculate  natural  ventilation 
pressures  differed  little  from  static  models  of  the  chimney 
effect  in  caves.  Beard  (1920)  employed  the  formula 

P  =  (1.3273)  (B)  (H)  ((1/R[ext  ])  -  ( 1/R[  int  ]) ) 

where  P  is  the  pressure  difference  (lbs/ft2) ,  B  is  the 
atmospheric  pressure  (in.  of  Hg) ,  H  the  height  of  the 
columns  (f t) ,  and  R£ext]  and  R[int]  the  average  temperature 
(°R)  of  the  external  and  internal  columns.  This  is 
equivalent  to  the  static  model  of  Trombe  (see  page  8) ,  with 
the  exception  that  temperature  rather  than  density  is 
averaged  over  height.  Penman  and  Penman  (1947)  adopt  the 
same  approach  although  their  pressure  differential  is 
couched  in  terms  of  a  motive  column.  This  motive  column  may 
be  visualised  by  assuming  that  the  two  static  columns  have 
the  same  density,  but  unequal  heights.  The  pressure 
differential  between  the  two  columns,  then,  is  merely  this 
common  density  times  the  difference  in  height  of  the  two 
columns,  times  the  acceleration  due  to  gravity. 


The  natural  ventilation  pressure  refers  only  to  the 
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static  pressure  difference  between  two  columns  of  air. 
Should  air  circulation  be  allowed  between  the  two  columns 
(that  is,  between  the  mine  or  cave  and  the  surface) ,  the 
ensuing  air  velocity  may  be  related  to  the  natural 
ventilation  pressure  by  the  following  formula,  which  is 
ascribed  to  Atkinson  (1886,  as  cited  in  Weeks,  1926) : 

(P)  (A)  =  (K)  (S)  (V)  2  , 

where 

P  =  pressure  difference  in  lb/ft2 
A  =  cr css-sectional  area  of  passage  in  ft2 
S  =  rubbing  surface  in  ft2 
V  =  air  velocity  in  ft/min 
K  =  coefficient  of  friction 

The  value  of  K  is  dependent  upon  the  units  used  and  includes 
the  average  density  of  the  air.  The  equivalent  formula  using 
metric  units  (i.e.,  the  appropriate  friction  coefficients) 
may  be  found  in  Stoces  (1954)  • 

This  static  approach  to  air  flow  in  mines  is  prevalent 
even  today.  McElroy  (1950)  and  Harris  and  Kingerly  (1973) 
make  minor  allowances  for  flow-induced  temperature  changes, 
but  their  models  are  by  no  means  dynamic.  McElroy 
acknowledges  the  effects  of  moisture  on  air  density,  but 
considers  its  inclusion  into  his  model  unwarranted.  Harris 
and  Kingerly  give  the  following  formula  for  the  density  of 


moist  air: 
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Density  =  (1 . 327)  ( 1/T)  (B-0. 378£) 

where  the  density  is  in  lb/ft3;  T  is  the  temperature  in 
degrees  Rankine;  B  the  atmospheric  pressure  in  inches  of 
mercury;  and  f  the  vapor  pressure  (in  inches  of  mercury)  at 
the  dewpoint  temperature.  They  then  state  that  temperature 
and  pressure  are  too  variable  for  this  sort  of  accuracy,  so 
that  the  vapor  pressure  can  be  ignored.  Their  model  then 
becomes  identical  to  that  of  Beard  (1920) ,  as  previously 
mentioned . 

Not  all  mine  engineers  follow  this  school  of  thinking, 
however.  The  departure  from  static  modelling  began  with  F.B. 
Hinsley  (1938,  as  cited  in  Williams,  1960;  1943a, b,c,d) . 
Hinsley  recognised  natural  ventilation  as  a  problem 
involving  the  thermodynamics  of  compressible  fluid  flow. 
Heat  added  to  (taken  frcm)  the  mine  air  is  balanced  by  the 
work  done  by  (on)  the  air  and  by  the  change  in  internal 
energy  of  the  air.  The  work  done  includes  changes  in  kinetic 
and  potential  energies,  changes  in  pressure  energy 
(Hinsley's  terminology) ,  and  the  work  done  against  friction. 
The  heat  added  includes  the  net  external  heat  added,  the 
frictional  heat,  and  any  mechanical  energy  added  (all 
variables  are  measured  per  unit  weight) .  Thus  the  governing 
energy  eguation,  between  any  two  cross* sections  i  and  j,  is: 

H  (1,2)  ♦  M  (  1 , 2)  ♦  F  (1,2) 

=  [Z(2)-Z(1)  ]  ♦  [{U(2)2  -  U(1)2}/2G] 

♦  [  P  (2)  V  (2)  *  P(1)V(1)  ]  ♦  F(1,2) 


, 
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+  £  E  (2)  -E  (1)  ] 


where 


H(i,  j) 

= 

net  heat  added  between  i  and  j 

M(irj) 

= 

mechanical  energy  added  between  i  and  j 

F(i,j) 

— 

frictional  heat  added  between  i  and  j 

I(  i.j) 

= 

work  done  against  friction  between  i  and  j 

Z(i) 

= 

potential  energy  at  i 

U(i) 

— 

air  velocity  at  i 

P(i) 

— 

pressure  at  i 

V(i) 

= 

specific  volume  at  i 

E(i) 

— 

internal  energy  at  i 

G 

— 

acceleration  of  gravity 

(Care  should  be  taken  to  assure  compatibility  of  dimensions 
when  mixing  heat  and  mechanical  energy  terms,)  In  the  case 
of  natural  ventilation,  there  is  no  mechanical  energy 
involved  (H(1,2)  =  0)  and  the  work  done  against  friction  is 
assumed  to  be  exactly  the  heat  gained  from  friction  (that 
is,  there  is  no  loss  of  heat  to  external  sinks;  F(1,2)  = 
F  (1,2))  •  If 

tf(i,j)  =  work  done  on  or  by  the  fluid  between  i  and  j  , 

the  preceding  equation,  ignoring  variations  in  potential  and 
kinetic  energy  reduces  to 

H(1,2)  +  F  (1,2)  =  W  (1 ,2)  ♦  [E(2)-E(1)  ]. 


Hinsley  showed  that  the  work  done  was 


' 


,  •  i  '  ' 

■ 


; 


l  '  (  I  <  ,  • 


18 


"2 

W(1,2)  =  PdV  , 

so  it  is  necessary  to  know  the  variation  of  pressure  with 
specific  volume  to  be  able  to  calculate  this  work. 

If  an  ideal  mine  is  assumed  (as  pictured  in  Figure  5) , 
then  air  moving  through  the  mine  (for  whatever  reason)  goes 
through  the  following  thermodynamic  cycle: 

1)  Air  enters  the  mine  and  descends,  gaining  heat  at 
the  expense  of  potential  energy. 

2)  At  the  bottom  of  the  downcast  shaft,  air  passes 
through  the  workings  and  gains  heat  from  the  strata. 

3)  The  air  rises  up  the  upcast  shaft,  losing  internal 
energy  to  potential  energy. 

4)  The  air  at  the  top  of  the  upcast  shaft  is  returned 

to  the  state  of  the  air  at  the  top  of  the  downcast 

shaft . 

On  a  pressure-specific  volume  (P-V)  graph,  this  cycle  is 
ideally  that  of  Figure  5.  The  net  work  done  per  lb  of  air 
(the  area  enclosed  by  the  cycle) ,  times  the  average  value  of 
the  air  density,  yields  the  natural  ventilation  pressure. 

Several  ether  authors  have  adopted  the  thermodynamic 
approach;  among  them,  Williams  (1960) ,  and  Hall  (1953, 

1967) .  While  the  approach  is  basically  correct,  there  are 
numerous  problems  associated  with  the  actual  calculations. 
If  actual  field  measurements  are  taken  and  plotted  on  a  P-V 
diagram,  the  determination  of  the  natural  ventilation 

pressure  is  reasonably  straightforward.  To  calculate  the  NVP 
from  surface  conditions  alone,  without  actually  measuring 
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the  pressure-volume  variation,  is  considerably  more 
difficult.  Hinsley  (1 943a, b, c,d)  assumed  that  flow  in  the 
shafts  was  frictional  adiabatic  (Hinsley* s  terminology  - 
adiabatic  modified  by  frictional  effects) ,  and  that  the 
collars  of  the  upcast  and  downcast  shafts  were  at  the  same 
elevation  sc  potential  energy  differences  would  sum  to  zero. 
Hall  (1967)  looks  extensively  at  these  assumptions  and 
attempts  to  analyse  the  effects  of  removing  them.  Hall  shows 
that  the  ideal  mine  cycle,  with  no  friction,  with  collars  at 
the  same  elevation,  with  no  moisture,  and  with  no  net  loss 
of  kinetic  energy  through  the  cycle,  is  actually  the  Joule 
or  Brayton  cycle.  This  cycle,  which  consists  of  two 
reversible  adiabatic  processes  joined  by  two  isobaric 
processes,  is  pictured  in  Figure  5  on  a  P-V  diagram.  If 
these  factors  are  taken  into  consideration,  the  cycle  should 
look  like  Figure  6,  but  Hall  is  unable  to  evaluate  their 
simultaneous  effects  and  must  resort  to  a  step-by-step 
evaluation  resulting  in  a  cycle  like  that  of  Figure  7. 
Assumptions  are  also  made  in  Hinsley's  model  concerning  the 
character  of  the  heat  transferral  in  both  the  shafts  and  the 
workings,  and  Hinsley  says  very  little  concerning  mine 
entrances  at  different  elevations,  the  situation  most 
analogous  to  dynamic  caves.  Clearly  improvements  can  be  made 
to  both  Hinsley *s  and  Hall's  models. 
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IV  -  THE  CCflPEESSIBlEMCDjl 


The  Approach 

Much  of  the  difficulty  with  mine  ventilation  studies 
involves  the  necessity  for  calculating  the  natural 
ventilation  pressure  despite  the  inability  to  integrate  all 
the  involved  factors  into  a  continuous  process.  While  many 
spelecmicroclimatologists  have  followed  the  approach  of 
calculating  the  NVP  and  then  the  velocity  of  the  air  flow 
from  the  Atkinson  Formula,  this  is  circular  reasoning, 
because  the  velocity  of  the  air  determines  largely  the 
pressure  variation  that  exists  (due  primarily  to  the 
variation  of  temperature  with  velocity  as  will  be  seen 
later) .  A  more  valid  approach  would  aim  at  calculating 
velocities  first.  For  the  speleologist,  this  would  be 
enough;  the  cave  atmosphere  is  determined.  For  the  mine 
engineer,  once  the  mine  portion  of  the  P-V  cycle  is  known, 
calculation  of  the  surface  parameters  completes  the  cycle 
and  the  NVP  may  then  be  calculated. 

The  major  obstacle  to  integrating  all  the  factors 
involved  in  the  calculation  of  air  velocities  is  the  problem 
of  estimating  the  net  heat  exchange  between  the  air  and  the 
rock  strata.  However,  as  previously  pointed  out,  this  has 
been  solved  theoretically  by  Wigley  and  Brown  (1971).  It  is 
therefore  possible  to  determine  the  parameters  of  the  cave 
atmosphere  under  steady-state  conditions,  given  the  external 
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atmospheric  conditions  at  the  two  entrances  and  certain  non¬ 
dynamic  properties  of  the  particular  cave.  The  approach  is 
via  the  thermodynamics  of  compressible  fluid  flow  (or  gas 
dynamics) ,  and  the  basic  model  is  described  in  the  next 
section. 

The  Model 

To  apply  the  techniques  of  gas  dynamics  to  an  actual 
cave  would  dc  little  to  elucidate  the  principles  involved. 
Variations  in  passage  gradient  and  cross-sectional  area,  and 
the  effects  of  passage  bending  would  serve  only  to  mask  the 
actual  chimney  process.  To  avoid  such  difficulties,  the 
passage  will  be  assumed  to  have  a  constant  radius  and 
gradient;  but  in  order  to  give  the  reader  some  insight  as  to 
their  influence,  these  assumptions  will  not  be  immediately 
introduced.  The  effects  of  bends  and  of  momentum  losses  due 
to  rapid  changes  in  passage  size,  direction,  or 
configuration  will  at  no  time  be  considered.  For  information 
in  this  area,  the  reader  is  referred  to  Shapiro  (1953),  from 
whom  much  of  the  following  is  taken,  fiadial  variations  in 
the  parameters  are  of  no  importance  here,  so  the  flow  will 
be  treated  as  one-dimensional  with  all  parameters  being 
averages  (for  a  more  precise  definition  of  average,  see 
Wigley  and  Brown,  1971)  over  each  cross-section.  For  the 
present,  only  the  steady-state  case  will  be  investigated. 

The  effects  of  moisture  in  the  air  deserve  special 
mention.  In  most  situations,  moisture  will  be  added  to  (or 
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taken  away)  from  the  air  by  condensation  (or  evaporation) • 
This  creates  a  substantial  heat  source  (or  sink)  which 
cannot  be  neglected.  Density  differences  due  to  moisture, 
however,  are  minor  and  will  be  ignored,  so  that  the  mass 
rate  of  flow  is  constant. 

The  Dynamics  and  Thermodynamics 

The  ideal  dynamic  cave,  with  constant  gradient  and 
cross-section,  is  pictured  in  Figure  8,  and  all  relevant 
parameters  are  explained  in  Table  1 .  The  basic  procedure  is 
to  look  at  an  incremental  section  of  passage  (called  the 
control  volume  -  see  Figure  9)  of  length  dx,  and  describe 
the  change  in  density,  pressure,  and  air  velocity  in  terms 
of  the  temperature  and  the  friction  along  the  walls.  The 
solution  is  obtained  from  the  simultaneous  solution  of  the 
equations  of  state,  continuity,  and  momentum.  For 
simplicity,  the  concepts  of  Mach  number  and  stagnation  state 
are  introduced;  and,  for  completeness,  the  energy  equation, 
although  not  required  for  the  solution,  is  also  derived. 

For  a  perfect  gas,  which  dry  air  approximates 
(corrections  for  moist  air,  to  avoid  a  variable  gas 
constant,  will  be  made  later) ,  the  equation  of  state  is 

P  =  E  (Eg)  T  ,  (4-1) 


from  which  logarithmic  differentiation  yields 
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dP  =  dD  ♦  dT 

P  D  T  (4.2) 

The  Mach  number  is  defined  to  be 

r  1*0.5 

M  =  V  |k  (Eg)  T  | 

L  J  (4.3) 

Thus  by  logarithmic  differentiation, 
dM£  =  dV2  -  dT 

M2  V2  T  (h.4) 

since  k  is  constant  for  a  perfect  gas.  The  mass  rate  of  flow 
is 

w  =  DAV  ,  (4.5) 

so  again  by  logarithmic  differentiation, 

0  =  dD  +  dA  +  dV 

DAV  (4.6) 

since  w  is  constant. 

The  momentum  equation  is  somewhat  more  complicated. 

"The  net  force  acting  on  the  material  within  the  control 
surface  is  equal  to  the  increase  of  momentum  flux  of  the 

streams  flowing  through  the  control  surface"  (Shapiro,  1953, 
p.  224) .  Thus,  assuming  the  divergence  of  the  walls  is  small 
(angle  of  divergence  approximated  as  0)  ,  we  have 

PA  +  PdA  -  (P  ♦  dP)  (A  +  dA)  -  t£  wall  ]dA[  wall  ]  -  dX 


w(V  ♦  dV)  -  w  (V) 


(4.7) 
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The  wall  friction  factor  is  defined  to  be 


f 


2tf  wall  1 

DV2 


w 


(4-8) 


so 


t[wall]  =  f  DV2 

2  (4. 9) 

Further,  by  defining  B  as  the  mean  hydraulic  diameter  (the 
cross-sectional  area  of  the  passage  divided  by  the  perimeter 
of  the  cross-section) ,  we  have 

£A£wallJ  =  4dx 

A  B  (4.10) 

Thus ,  from  (4.7),  (4.9),  and  (4.10),  and  neglecting  second 

order  terms, 

wdV  =  -AdP  -  2ADV 2f dx  -  dX  . 


B 


(4.11) 


But,  from  (4.1)  and  (4.3), 


DV2 


kPM2  , 


(4.12) 


so  from  (4.5),  (4.11)  and  (4.12),  we  get 


kPAM2dV 

V 


AdP  -  2 k P A M£f dx  -  dX 
B 


(4.13) 


Dividing  by  PA  and  rearranging. 


-kM2dV 

V 


k M  2 2 f d x  -  dX 
B  PA 


dP 

P 


(4. 14) 
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But  X  is  the  component  of  weight  working  against  the  flow, 
so  if  Z[ x  ]  is  the  height  above  the  lower  entrance  at  x 
distance  into  the  cave, 

dX  =  BAGZ • [ x  ]dx  ,  (4.15) 

and  thus,  frcm  (4.14),  (4.15),  and  (4.1), 

dP  =  -kM2d V  -  kM2 2f dx  -  DGZ'fx  Idx 
P  V  B  P 

=  -kM2dV  -  kM2 2fdx  -  GZl£xJdx 

V  B  (Eg)  T  (4.16) 

The  energy  equation,  while  not  required  for  solution  of 
the  problem,  is  nevertheless  instructive.  For  flow  through 
the  control  surface,  the  rate  of  heat  flow  is 

r  i 

wdQ  =  w  (h  +  dh)  -  wh  +  j  V£  +  dV2  Jw 

I  2  2  1 

L  J 

-  wV£  -  wc[p]I(z  +  dz) 

2 

♦  wc[p]Lz 

=  wdh  ♦  wdV2  -  wc[p]ldz  , 

2  (4.17) 

where  dQ  is  the  net  heat  per  unit  mass  added  to  flow  by 
conduction  cr  radiation  from  external  sources.  Dividing  by 
w,  we  get 

dQ  =  dh  dV2  -  c[p]Ldz 

2  (4.18) 


We  may  break  dh  down  to 
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dh  =  -dh[  pr  ]  +  c[  p  ]dT  ,  (4.19) 

yielding,  upon  rearrangement, 

d Q  +  dh[  pr  ]  +  c£p  ]Ldz 

=  c[  p  ]dT  ♦  dyz 

”2  (4.20) 

Division  by  c£p]T  gives 

d£_  +  _ dh£ pr J  ♦  cf p  )Ldz  =  dT  + dV£_  , 

c£p]T  T  2c£p]T  (4.21) 

and  substitution  from  (4.4)  and  the  relationship  between  Eg, 
c[p],  and  k  yields 

dQ  +  dhfprl  +  cf  plLdz  =  dT  +  (k-1 )  M2  dV£ 

c£p]T  T  2  V2  (4.22) 

Thus,  the  heat  input  from  the  strata  plus  the  heat  input  (or 
output)  from  condensation  (evaporation) ,  plus  the  energy 
lost  to  the  gain  in  elevation  (negative) ,  is  described 
totally  by  the  temperature  and  velocity,  as  would  be  be 
expected . 

It  is  advantageous  here  to  introduce  the  concept  of 
stagnation  state,  which  is  defined  as  that  state  arrived  at 
by  isentropically  reducing  the  flow  to  zero  velocity.  From 
(4.18) ,  neglecting  potential  energy  and  phase  changes, 

ho  =  h  +  V£  , 

2  (4.23) 

from  which  the  following  relationship  can  be  derived: 
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To  =  T  ♦ _ y  2 

2c£p] 

From  (4,3) ,  this  may  be  written 


(4.24) 


r  t 

To  =  T  |  1  +  (k-  1)  M2  | 

I  2  | 

L  J 


Similarly, 


Po 


r  Tk/(k-1) 

P I  1  ♦  JkrJIMi  j 

I  2  | 

L  J 


Do 


r  tV(Jc-I) 

E  |  1  +  (k- 1 )  Mg  J 

t  2  j 

L  J 


(4.25) 


(4.26) 


(4.27) 


fiather  than  solve  directly  for  density,  pressure,  and 
velocity,  it  is  easier  to  solve  for  the  Mach  number 
variation  along  the  passage  (to  be  discussed  later) .  Then, 
from  equations  (4.2),  (4.3),  (4.4),  and  (4.6)  and  the 

variation  of  the  temperature  along  the  passage,  the  density, 
pressure,  and  velocity  variations  and  their  related 
stagnation  properties  may  be  obtained. 


The  solution  for  density  from  (4.4)  and  (4.6),  assuming 
the  Mach  number  and  temperature  variation  to  be  known,  is 


dD 

D 


dA  -  dV 
A  V 

dA  -  dVf 
A  2V2 


v  * 

4 


* 


■  ' 


' 

*  (  t 


» 


' 

»  , 

. 


28 


=  -  d  A  -  dM£  -  dT 

A  2M2  2T 


(4.28) 


If  (4.28)  is  integrated  between  any  two  sections  1  and  2  and 
exponentials  taken,  we  arrive  at 


r 

t  0.5 

D(2) 

=  iJll  MJJl  1 

Till  1 

D(1> 

4(2)  M(2)  | 

T  (2)  | 

L 

J 

(4.29) 

Similarly 

,  for  pressure. 

from  (4.2)  we  have 

PJ2L 

=  SJ2±  TJ21 

P(1) 

E(1)  T  (1) 

• 

(4.30) 

But  substitution  from  (4 

.29)  gives 

r 

i  0 .5 

PJ21 

=  Min.  i 

TJ21  | 

P(1) 

4(2)  »(2)  | 

T  (1)  1 

L 

J 

(4.31) 

The 

solution  for  Mach  number  variation  is 

somewhat  more 

complicated.  Starting 

with  equations  (4.4), 

(4.6)  ,  and 

(4.2),  we 

get 

dM2 

=  d  V2  -  dT 

M2 

V2  T 

=  -2dD  -  2dA  - 

dT 

D  A 

T 

=  -2dP  -  2d A  + 

dT  . 

P  A 

T 

Substitution  from  the  momentum  equation,  (4.16) 

,  gives 

d  M£ 
M2 


kM2dV 2  +  kM 2 4fd  x  ♦  2GZT  x  Idx  -  2dA  +  dT 
V2  B  (Eg) T  AT 


< 


«  -  _  ’ 


1  _ 
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and  from  (4.4), 


M£  =  k«24M£  +  (1  «■  kM2 )  dT  +  kM24fdx 

M2  M2  T  g 

♦  2GZ^£xJdx  -  2d A  . 

(Rg)T  A 


By  collecting  dM2  ,  this  yields 

M2 


(1  -  kM2 ) dM£ 

~  M2 

=  (1  +  kM2) dT  +  kM24fdx  *  2GZ'f  x  ]dx  -  2dA  . 

T  B  (Bg)  T  A 

r  i  r  t 

dM2  =  |  1  +  kM2  |  dT  ♦  |  _4kM£ |  fdx 

M2  |  1  -  kM2  |T  |  1  -  kfl2  |  B‘ 

L  i  l  J 


♦ 


r 


i 


I  _ 2 _ | GZ«( x  ]dx 

|  1  -  kM2  |  (Eg) T 

L  J 


r  i 

| _ 2 _ |  dA 

|  1  -  kM2  I  A 


L 


J 


(4.32) 


Only  under  very  specific  conditions  may  this  equation  be 
solved  analytically.  For  a  numerical  solution,  we  first 
invoke  the  simplif ications  previously  mentioned.  The  cross- 
sectional  area  becomes  constant  with  constant  radius  E 
(necessitated  by  the  temperature  solution  discussed  later) . 
Thus  dA/A  =  0  and  B  =  2E.  Further,  the  gradient  is  constant, 
so  Z[ x  ]  =  xsinG  and  Z*[x]  =  sinG.  With  these 
simplifications,  equation  (4.32)  becomes 
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r  i  r  i 

=  I  1  dT  +  |  2kflz _ j  fdx 

M2  |  1  -  kM2  |  T  |  1  -  kH2  |  B 

t  J  L  . 


r  i 

*  I  _ 2 _ j  GsinOdx 

I  1  -  kM2  |  (fig) T 
*-  J 


For  manipulative  simplicity,  let 


FI  =  M2(1  +  k M2)  , 

1  -  kM2 


FII  =  _2kM* _  , 

1  -  kM2 


Fill  = _ 2M2 _  . 

1  -  kM2 


Then  (4,33)  becomes 


dM2  =  FIdT  ♦  FIlfdx  +  FIIIGsinOdx  . 

T  B  ~(Rg7  T 


(4.33) 


(4.34) 

(4.35) 

(4.36) 


(4-37) 


By  assuming  that  the  coefficients  of  dT  and  dx  are  constant 
at  their  mean  value  over  the  domain  of  integration 
(trapezoidal  rule) ,  equation  (4.37)  may  be  integrated 
between  sections  1  and  2  to  yield 


\  11121  ♦  2I-U1  ]  I  T(2)  -  T(1) 
I  T(2)  T(1)  ||  2 


L. 


J  L 


T 

I 

I 

J 


M (2)  2  -  M  (1) 2 


.  ) 
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I  I  Lilli  2] _ ±_l  IIJ11  j 

I  R 


1  r 

I  1 
I  I 

J  L 


x!2l_r. 

2 


.xill  I 
I 


I  Fill  (2)  +  FI II J 11  (Gsinf 

I  1(2)  T(1)  |  (fig) 


I  x(2)  -  x(1)  | 


I 


I 


(4.38) 


Now  from  equations  (4.25),  (4.26),  (4.27),  (4.30),  (4.31), 
and  (4.38),  a  complete  solution  for  density,  velocity, 
pressure,  and  the  related  stagnation  properties  is  available 
for  any  cross-section  (2)  in  terms  of  some  previous  cross- 
section  (1)  and  the  temperature  distribution.  As  FI,  FII, 
and  Fill  all  contain  M(2),  the  solution  for  M (2)  in  equation 
(4.38)  will  of  necessity  be  an  iterative  one. 


Temperature  Variation  -  Dry  Air 


For  illustrative  purposes,  let  us  initially  look  at  the 
case  where  there  is  no  moisture  whatsoever  in  the  cave  or 
free  atmosphere.  In  this  instance,  the  variation  of 
temperature,  assuming  fully— developed  turbulent  flow,  is 
expressed  (in  the  case  of  a  horizontal  passage)  by  the 
differential  equation 


dT  =  JTHALL_2_I1 
dx  xo 


(4.39) 


(see  Shapiro, 


1953;  Wigley  and  Brown,  1971).  Wigley  and 
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Brown  (1971)  have  shown  that,  for  smooth-walled  passages, 
xo  =  IJPrlJRei  , 

2  "(Nu)  (4.40) 

where  Pr  is  the  Prandtl  number.  Re  the  Reynolds  number,  and 
Mu  the  Nusselt  number  of  the  flow.  Empirical  relationships 
have  been  worked  out  (see  for  instance  Kays,  1966)  which 
yield  the  Nusselt  number  in  terms  of  the  Prandtl  and 
Reynolds  numbers,  but  the  correct  relationship  depends  upon 
the  nature  of  the  heat  flew  between  rock  wall  and  air.  It  is 
necessary,  therefore,  to  look  at  what  determines  the 
temperature  of  the  rock  wall. 

To  analyse  TWALL,  consider  what  conditions  would  exist 
within  the  solid  rock  mass  in  the  absence  of  a  cave.  "The 
thermal  condition  of  a  rock  mass  at  a  given  time  is  a 
function  of  its  conductivity,  density,  heat  capacity,  rate 
of  internal  heat  generation  and  loss,  temperature 
fluctuations  of  the  surface,  and  the  regional  geothermal 
gradient"  (Lange,  1954,  p.21) .  By  assuming  a  level-surfaced 
homogeneous  rock  mass  of  constant  thermal  conductivity  and 
diffusivity,  Lange  has  calculated  the  variation  of 
temperature  with  depth  from  a  sinusoidal  variation  in 
surface  temperature.  If  k  is  the  thermal  conductivity  of  the 
rock,  and  d  and  c[ p  ]  its  density  and  heat  capacity  (constant 
pressure) ,  then  the  temperature  6  at  depth  z  and  time  t 
resulting  from  a  surface  temperature  fluctuation  of  Ncoswt 
is 
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G(z,t)  =  Nexp  (-mz)  cos  (wt  -  mz)  , 


(4.41) 


where 


m 


r  -i0.5 

I  wdcf  p  1  | 

I  2k  | 

i  j 


(4.42) 


This  fluctuation  is  superimposed  upon  a  mean  geothermal 
gradient  which  is  assumed  to  be  linear  with  depth.  The 
period  of  the  surface  fluctuation  is  2  tt/m  and  the  amplitude 
is  N,  and  the  thermal  diffusivity  k/dc[p]  depends  upon  the 
particular  rock.  For  certain  marbles  that  Lange 
investigated,  the  thermal  diffusivity  was  0.00502  m2/hr,  so 
using  this  value,  a  diurnal  temperature  variation  of  20°C 
would  yield  a  subsurface  variation  of 

G(z,t)  =  lOexp  (-5. 1  z)  cos  (.262t  -  5.1z)  .  (4.43) 

From  equation  (4.43),  by  setting  (.262t  -  5.1z)  =  0,  the 
depth  at  which  a  2°C  variation  exists  (as  a  result  of  the 
20°C  surface  variation)  is  0.45  meters.  Clearly,  daily 
surface  variations  will  have  no  significant  effect  upon  rock 
temperatures  at  depth. 


Now  let  us  assume  a  40°C  annual  variation.  Equation 
(4.41)  now  becomes 


0(z,t)  =  20exp (-0.267z) cos ( .0007 17t  -  0.267z)  , 

(4.44) 


so  the  depth  at  which  a  2°C  variation  is  felt  is  11.2 
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meters.  Again,  seasonal  effects  are  minimal  with  depth.  This 
process  may  be  continued  indefinitely,  and  it  becomes 
reasonably  cbvious  that  only  cycles  of  very  long  duration 
can  possibly  affect  rock  temperatures  to  any  depth.  Thus, 
for  our  purposes,  the  geothermal  gradient  may  be  considered 
time-independent. 

Despite  the  lack  of  dependence  on  temporal  variations, 
the  determination  of  rock  temperature  is  still  complex.  In 
low-lying  karst  areas,  such  as  central  Kentucky,  the  problem 
is  not  too  difficult;  the  temperature  at  depth  may  be 
extrapolated  from  the  average  surface  temperature  and  the 
regional  geothermal  gradient,  the  latter  usually  being 
between  0.01  and  0.05°C/meter  (see  Carslaw  and  Jaeger,  1959; 
as  cited  in  Ford  et  al. ,  1975) .  In  alpine  areas,  however, 
the  change  of  mean  annual  temperature  with  elevation,  and 
especially  the  effects  of  a  highly  irregular  surface, 
severely  complicate  determination  of  the  rock  temperature  at 
depth,  and  such  an  endeavor  would  be  beyond  the  scope  of 
this  thesis.  One  approximation  would  be  to  assume  the 
surface  was  level  and  extrapolate  from  the  local  mean  annual 
temperature  and  the  regional  geothermal  gradient  (see  Figure 
10)  .  The  error  should  not  be  large  if  the  local  surface  is 
not  steep  and  the  relief  not  great.  This  approach  is  common 
in  mine  engineering,  where  surface  irregularities  are  often 
small  in  ccmparison  to  mine  depths.  Caves,  however,  are 
considerably  mere  restricted  in  their  location;  that  is, 
they  are  formed  by  groundwater  flow  and  hence  do  not  extend 
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significantly  deeper  than  the  lowest  point  at  which  the 
limestone  outcrops  (at  the  time  of  their  formation) ,  which 
is  the  approximate  point  at  which  their  waters  resurface. 
Thus,  cave  depths  are  generally  of  the  same  order  of 
magnitude  as  (or  less  than)  the  local  topographic  relief, 
and  the  variation  of  mean  annual  surface  temperature  with 
elevation  significantly  affects  the  cave  rock  temperature  at 
depth.  A  possible  alternative  approximation  would  be  to 
assume  that  the  temperature  of  the  rock  varies  linearly  with 
the  distance  normal  to  the  mean  surface  (see  figure  10) . 
Many  caves  descend  or  ascend  at  approximately  the  same 
gradient  as  the  surface;  hence  estimates  of  the  rock  wall 
temperature  by  this  method  would  yield  reasonable  results. 

The  simplest  approximation  is  to  assume  an  isothermal 
rock  mass.  For  shallower  caves  in  low-lying  areas,  this  is  a 
reasonable  approximation.  For  caves  in  alpine  areas  or  at 
great  depth  in  low-lying  regions,  however,  the  isothermal 
case  can  only  be  considered  a  first  approximation.  Those 
interested  in  further  approximation  methods  are  referred  to 
Kays  (1966),  where  an  in-depth  discussion  of  the 
relationship  between  the  Nusselt,  Prandtl,  and  Reynolds 
numbers  (for  flow  through  constant  diameter  pipes)  is  found. 
For  the  purposes  of  this  paper,  only  the  isothermal  case 
will  be  investigated,  because  of  the  mathematical  simplicity 
it  affords. 


The  foregoing  discussion  concerned  rock  temperatures  at 
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depth  in  a  homogeneous  rock  mass.  The  introduction  of  a 
dynamic  cave  into  this  rock  mass  necessarily  alters  somewhat 
the  temperature  distribution  within  the  rock,  for  not  only 
does  the  rock-wall  temperature  affect  the  internal  air 
temperature,  but  also  the  reverse  is  true.  The  effects  are 
most  noticeable  in  the  entrance  area,  where  the  internal  air 
temperature  differs  most  from  the  temperature  of  the  rock 
mass.  Beyond  the  entrance  zone,  however,  the  effects,  should 
be  minimal  because  the  temperature  difference  rapidly 
becomes  quite  small. 

For  the  sake  of  mathematical  tractability,  Wigley  and 
Brown  (1971)  ignore  heat  transfer  within  the  wall  and 
replace  it  with  the  boundary  condition  that  the  wall 
temperature  TWALL  is  constant.  Under  these  somewhat 
stringent  conditions,  the  relationship  between  the  Nusselt, 
Prandtl,  and  Beynolds  numbers  (for  a  constant  pipe  diameter) 
is 


0.6  0.8 

(0.021)  (Pr)  (Be) 


(4.45) 


Nu 


(see  Kays,  1966;  as  cited  in  Wigley  and  Brown,  1971) .  This, 
with  equation  (4.40),  yields 


xo 


0.4  0.2 

(24)  (Pr)  (Re)  B 


(4.46) 


From  the  definitions  of  Pr  and  Be,  this  in  turn  gives 
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r  t  0.  2 

xo  =  (28)  |  (coeff  of  viscosity)  cf  p  I2  | 

|  (thermal  conductivity) 2  [ 

L  J 

1.2  0.2 

•  (R)  (DV)  .  (4.47) 

While  D  and  V  vary  along  the  passage,  DV  is  constant  and 
equal  to  w/A.  Further,  for  dry  air, 

c[p]  =  1004  m2/sec2-°C  (Hess,  1959), 

and  since  the  coefficient  of  viscosity  and  the  thermal 
conductivity  vary  little  with  temperature  and  density,  a 
reasonable  approximation  is  gained  by  considering  them 
constant.  At  10°C,  their  values  (by  linear  interpolation 
from  Weast,  1970)  are: 

-5 

coefficient  of  viscosity  =  1.774(10)  kg/m-sec  , 

-2 

thermal  conductivity  =  2.482(10)  J/m-sec-°C 

Thus,  equation  (4.47),  considering  that  A  =  ttr2  ,  becomes 

0.8  0.2 

xo  =  (170)  (R)  (w)  .  (4.48) 

The  previous  formula  assumes  the  passage  is  smooth-walled. 
If  this  is  net  the  case,  xo  should  be  corrected  with  use  of 
the  actual  friction  factor  (see  Kays,  1966),  as  follows: 
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r  t0.5  0.8  0.2 

xo  =  I  i Csmooth  1  |  (169.86)  (R)  (w) 

I  f  I 

1  J  (4.49) 

Since  TWALL  is  constant  for  the  isothermal  rock  mass 
considered  in  this  paper,  and  xo  is  constant  for  steady  flow 
conditions,  (assuming  a  constant  diameter  pipe)  equation 
(4.39)  may  now  be  solved.  The  solution  (given  T=TEXTL  at 
x  =  0)  is: 

T  =  TWALL  ♦  (TEXT!  -  TWALL) exp (-x/xo)  ,  (4.50) 

which  describes  a  simple  exponential  decay  of  the  entrance 
air  temperature  to  the  rock-wall  temperature,  with  distance 
into  the  cave. 

To  alter  this  formula  for  our  non-horizontal  passage, 
as  Wigley  and  Brown  (1971)  state,  the  appropriate  adiabatic 
lapse  rate  (in  this  case,  the  dry  rate)  need  only  be 
superimposed  as  a  non-dynamic  effect.  This  cannot  be  the 
case,  however,  for  any  change  in  temperature  due  to  a  change 
in  elevation  must  affect  subsequent  temperatures  in  the  cave 
(by  equation  (4.39)).  The  proper  method  of  taking  these 
factors  into  account  is  to  replace  the  actual  temperature 
with  the  potential  temperature  (Tp) ,  which  is  the 
temperature  after  the  process  is  adiabatically  reduced  to 
isobaric  conditions.  The  actual  and  potential  temperatures 


are  related  as  follows 
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r  R/c[  p  ] 

Tp  =  T|  1000  | 

I  P  I 

t  j 


It  may  be  shewn  (from  the  second  law  of  thermodynamics)  that 
eguation  (4.39)  becomes 


TdlnTjD  =  (THAU  -  Tj  . 
dx  xo 

so  that 


dT  -  RTdlnP  =  (TWALL  -  . 

dx  c[ p  ]dx  ~  xo 


But 


-RTdlnP  =  _-RTdP_ 
c[p]dx  Pc[p]dx 

=  -GdP _ 

DGc[  p  ]dx 

=  LsinOdP 
~DGdz 


By  approximating  the  density  variation  with  height  as  that 
under  static  loads,  this  yields 

-RTdlnP  =  LsinO  , 
c[ p  ]dx 

so  we  see  that  dT/dx  in  the  level  (isobaric)  case  need  only 
be  replaced  by  (dT/dx)  +  IsinO  in  the  non-level  case. 
Solving  eguation  (4.39)  with  this  modification  yields 

T  =  TRAIL  -  LxosinO 

♦  (TEXTL  -  TWALL  -  LxosinO) exp  (-x/xo)  .  (4.51) 
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Temperature  Variation  -  Moist  Air 

The  introduction  of  moisture  to  the  air  complicates 
substantially  the  determination  of  temperature.  While  dry 
air  and  water  vapor  separately  approximate  perfect  gases, 
the  combination  does  not  (without  resort  to  a  variable  gas 
constant  dependent  upon  the  humidity  of  the  air) .  This  may 
be  overcome  by  introducing  the  concept  of  virtual 
temperature,  which  is  defined  as  "...the  temperature  at 
which  dry  air  would  have  to  be  in  order  to  have  the  same 
density  as  a  sample  of  moist  air,  assuming  both  have  the 
same  pressure"  (Hess,  1959,  p.60) .  The  use  of  virtual 
temperatures  allows  the  employment  of  the  eguation  of  state 
for  dry  air  under  all  circumstances.  The  general  formula  for 
virtual  temperature  is: 

Tv  =  1  (1  +  0 . 609g) 

Another  important  effect  related  to  the  introduction  of 
moisture  to  the  air  is  the  heat  addition  (or  loss)  from 
condensation  (or  evaporation) •  Wigley  and  Brown  (1971)  have 
also  taken  this  into  consideration  and  have  derived  a 
modified  form  of  the  differential  eguation  (4.39)  which 
accounts  for  condensation  and  evaporation.  This  modified 
equation  is: 
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r  t 

dT  =  J_|  TWALL  -  T  -  HE  (QWALL  -  QSPEC)  |  , 

dx  xo|  c[p]  | 

L  -•  (4.53) 


where 


dQSPEC  =  E (QWAL1  -  QS EEC)  , 

dx  xo  (4.54) 

and  E  is  the  fraction  of  the  passage  wall  that  is  moist. 
Generally,  E  would  vary  along  the  passage;  hence  it  would  be 
a  function  cf  x.  For  most  applications,  however,  it  is 
sufficient  tc  assume  that  E  =  1  Since  (4.53)  becomes 

(4.39)  when  E  =  0  ,  we  will  retain  E  in  the  formula  as  a 

constant,  sc  that  moist  and  dry  solutions  may  be  compared. 
The  solution  to  (4.54)  under  the  same  assumptions  as  for  dry 
air  (xo  and  TMLl  constant,  hence  QWALL  also  constant)  ,  is 

QSPEC  =  QWALL  +  (QEXT  -  QWALL)  exp  (-Ex/xo)  .  (4.55) 

Substituting  this  into  (4.53)  gives 


dT 

dx 


r  t 

1_ |  TWALL  -  T  +  HE  (QEXT  -  QWALL) exp (-Ex/xo)  |  , 

xo |  C[p]  I 

L  J 

(4.56) 


which  may  be  solved  to  yield 

T  =  TWALL  ♦  (TEXTL  -  TWALL)  exp  (-x/xo) 

♦  HEl^EXT_r_QKALLlxex_pl-Ex/xol  , 

c(  p  ]xo  (4.57) 


Again,  for  our  non- horizontal  passage,  this  becomes 
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T  =  THALL  -  Lxosin-0 

♦  (TEXTL  -  THALL  +  LxosinO) exp (-x/xo) 

♦  HE  (QEXT  -  QWALL) xexp  (~Ex/xo)  , 

c[p]xo  (4.58) 


Equation  (4.53)  does  not  always  apply,  however.  In 
certain  circumstances  supersaturation  would  be  predicted  by 
this  equation,  and  alterations  must  be  made  to  avoid  this. 
If  xsat  is  the  point  at  which  saturation  is  achieved,  and 
TSAT  the  temperature  at  that  point,  then  the  solution  (see 
Wrgley  and  Brown,  1971,  for  details)  is 

r  i 

(1  +  [  J  1  ])  ln|  TSAT  -  THALL  1  *  [  J 1  ][  J2  ]  (TSAT  -  T) 

|  T  -  THALL  | 

L  J 

=  x  -  xsat  , 

xo  (4.59) 


which,  rearranged,  becomes 


T 


=  THALL  ♦  (TSAT  -  THALL) 


r  i 

•  exp  |  -x  ♦  xsat  +  xof  J1  If  J2  1  (TSAT  ~  T)_  | 
|  (1  +  £  J1  ])  xo  | 

L  -« 


(4.60) 


where 


£J1]  = _ (QWALLLjHj _ , 

(Ev )  c£  p  ]  (THALL)  2 


and 
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[J2]  =  _ H _  . 

(Ev) (TWALI) 2 

Again,  equation  (4.60)  can  be  modified  for  the  non¬ 
horizontal  case  to  yield 

r  t 

T  =  TWALL  -  ac  +  |  TSAT  -  TWALL  ♦  _ac_  \ 

bc+1  |  bc  +  1  | 

L  J 

r  1 

•  exp  |  Jbcjljglxgat  -  x)  -  (T  -  TSAT)  (xo)b(bc+1)  |  , 

I  a  (xo)  | 

i  j 

(4.61) 

where 

a  =  1  +  [J1] 

b  =  £J1][J2] 

c  =1  (xo) sin6 

Thus  equation  (4.58)  is  used  when  x<xsat  and  equation  (4.61) 
when  x>xsat.  Both  equations  must  be  converted  to  virtual 
temperature  by  use  of  equations  (4.52)  and  (4.55)  if  they 
are  to  be  used  to  calculate  the  variation  in  Mach  number  in 
equation  (4.38).  For  x<xsat  the  formula  is 

Tv  =  T ( 1  ♦  0.609 (QWALL  ♦  (QEXT  -  QWALL)  exp  (-Ex/xo)  )  )  , 

(4.62) 


where  T  ccmes  from  equation  (4.58).  For  x>xsat. 

Tv  =  T  (1  +  0.609  (QSAT[T])  )  ,  (4.63) 

where  T  is  found  by  iteration  from  equation  (4.61). 
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The  description  of  the  temperature  distribution  is  now 
complete.  Equations  (4.58)  ,  (4.61),  (4.62),  and  (4.63)  may 
be  used  in  conjunction  with  equation  (4.38)  to  describe  the 
variation  of  Mach  number  along  the  passage  given  suitable 
boundary  conditions  to  be  discussed  later.  The  mass  rate  of 
flow  is  required  for  the  temperature  solution,  and  while  the 
stagnation  density  will  be  calculable  from  the  boundary 
conditions,  the  velocity  (and  hence  the  actual  density)  must 
be  found  iteratively.  The  friction  factor,  required  for  both 
the  description  of  the  temperature  variation  and  the  Mach 
number  variation,  is  discussed  below. 

The  Frictio n  Factor 

The  friction  coefficient  f  is  the  Fanning  friction 
factor,  which  relates  the  wall  shear  stress  to  the 
difference  between  the  actual  pressure  and  the  stagnation 
pressure  (commonly  called  the  velocity  pressure  -  see 
equation  (4.8)).  It  is  a  dimensionless  parameter  which  may 
be  shown  by  dimensional  analysis  to  depend  upon  the  Mach 
number,  the  Reynold's  number,  and  the  relative  roughness  of 
the  passage  wall  (Thompson,  1972).  The  dependence  of  f  on 
the  Mach  number,  however,  is  generally  small,  especially 
when  the  Mach  number  is  extremely  low  (<0.1),  as  is  the  case 
here.  Thus,  f  is  here  considered  to  be  a  function  only  of 
the  Reynolds  number  and  the  relative  roughness  (see 
Schlichting,  1968;  as  cited  in  Chapman  and  Walker,  1971)  . 
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for  the  steady  flow  case  (with  constant  pipe  radius)  being 
investigated,  the  Reynolds  number  is  treated  as  constant,  so 
that  f  depends  only  upon  the  relative  roughness.  The 
relative  roughness,  which  is  the  average  height  of  the  wall 
irregularities  divided  by  the  passage  diameter,  may  be 
expected  to  vary  along  the  passage.  Because  the  surface 
irregularities  are  only  average  values  over  a  distance,  a 
common  simplifying  assumption  is  that  the  friction  factor  is 
constant  for  a  given  problem  such  as  this  one.  Therefore, 
one  need  only  look  up  the  appropriate  value  of  the  friction 
factor  for  the  type  and  roughness  of  the  rock.  These  values 
are  commonly  found  in  any  modern  mine  engineering  handbook. 

In  hydraulics,  the  most  common  source  for  tables  of 
friction  factors  is  Moody  (1946) .  These  values,  however,  are 
for  the  Darcy  friction  factor,  rather  than  the  Fanning 
friction  factor.  The  Darcy  friction  factor  is  exactly  four 
times  as  large  as  the  Fanning,  and  the  reader  is  cautioned 
to  determine  which  a  given  set  of  tables  refers  to,  as  it 
often  is  not  stated. 

For  application  to  caves  and  mines,  however,  common 
engineering  sources  for  friction  factors  (such  as  the 
previously  mentioned  Moody  (1946))  are  generally  not 
sufficient,  because  the  range  of  relative  roughness  is  too 
low.  We  turn,  therefore,  to  the  mining  literature,  and 
should  immediately  caution  the  reader  as  to  what  will  be 
found.  A  multitude  of  definitions  exist  here  for  friction 
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factors.  They  are  all  denoted  by  K  and  seldom  can  one 
readily  determine  the  relationship  between  one  K  and 
another.  This  is  primarily  due  to  the  K's  being  dimensional 
while  f  is  dimensionless,  and  to  their  being  based  on  both 
the  Darcy  and  Fanning  fricticn  factors.  The  most  common 
friction  factor  in  mining  is  the  Atkinson  friction  factor, 
which  may  be  related  to  the  Fanning  friction  factor 
(assuming  the  average  weight  of  a  cubic  foot  of  air,  at  60°F 
and  atmospheric  pressure,  is  0.0764  lbf)  by 

f  =  3.03  K  , 


where 

the  K-tables 

have  been 

tabulated 

from 

velocity 

measurements  made  in 

thousands  of 

feet  per 

minute 

(Bunt, 

1960)  . 

Hartman  (1961) 

,  however,  relies  on  the 

relationship 

f 

=  1.214(10) 

7  K 

where 

the  K-tables 

have  been 

calculated 

from 

velocity 

measurements  in  feet  per  minute.  Conversion  to  thousands  of 
feet  per  minute  gives 

f  =  12.14  K  , 

from  which  we  see  that  Hartman's  friction  factor  K  must  be 
based  on  the  Darcy  friction  factor.  However,  he  makes  no 
distinction  between  Darcy  and  Fanning  friction  factors, 
which  seems  to  be  a  common  occurrence  in  mining  literature. 
The  reader  is  thus  again  cautioned  to  be  sure  of  the  meaning 
of  the  friction  factor  he/she  has  obtained. 
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In  determination  of  the  relaxation  length,  it  is 
necessary  tc  relate  the  actual  Fanning  friction  factor  of 
the  passage  to  the  friction  factor  of  the  same  passage  if 
the  walls  were  smooth  (relative  roughness  is  zero) .  This  has 
been  determined  empirically  as 

f[ smooth  ] 

r  r  n“2 

=  0. 25  |  0.08591n| _ _ Re _ |  j  (4.64) 

|  |  l79641n(Re)  -  3.8215  || 

i  i  j  j 

(Techo,  Tickner,  and  James,  1965;  as  cited  in  Chapman  and 
Kalker,  1971),  which  becomes,  under  the  previous  assumptions 
(page  37)  , 

f i smooth  ] 

r  r 

=  0.25|  0. 08591n| _ 35886. 12w/R _ I  j 

|  |  1 .9641n  (35886. 12w/R)  -  3.8215  j| 

l  L  JJ 

(4.65) 

Kith  the  temperature  variation  and  the  friction  factor 
determined,  the  model  is  now  complete  and  all  parameters  may 
be  determined. 

The  Procedure 

The  procedure  to  calculate  the  microclimatic  parameters 
reguires  multiple  iterations,  and  is  therefore  most  suited 
to  manipulation  using  a  digital  computer.  The  basic  premise 
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of  the  model  is  that  there  is  a  smooth  transition  through 
the  cave  from  the  atmospheric  pressure  at  the  lower  entrance 
to  the  atmospheric  pressure  at  the  upper  entrance.  It  is  the 
variation  of  air  velocity  that  maintains  this  transition.  It 
is  assumed  that  the  stagnation  pressure  and  temperature  at 
the  lower  entrance  are  the  atmospheric  pressure  and 
temperature  at  that  point,  since  nearly  all  the  energy  lost 
in  the  drop  of  pressure  and  temperature  as  the  air 
accelerates  into  the  cave  is  maintained  in  the  kinetic 
energy  of  the  air  (that  is,  the  process  is  isentropic)  . 
Further,  the  actual  pressure  at  the  upper  entrance  must 
balance  the  atmospheric  pressure  there.  Once  these  boundary 
conditions  are  set,  calculations  may  begin. 

1)  Specify  the  external  conditions:  atmospheric 

pressure  at  the  upper  and  lower  entrances  and 
temperature  and  relative  humidity  at  the  lower 
entrance  (for  winter  conditions) . 

2)  Specify  the  cave  passage  parameters:  length,  radius, 
angle  of  rise,  internal  wall  temperature,  and 
friction  coefficient. 

3)  Calculate  the  saturated  and  then  the  actual  specific 
humidity  for  atmospheric  conditions  at  the  lower 
entrance. 

4)  Estimate  the  initial  velocity  of  the  air  at  the 

lower  entrance. 

5)  Calculate  the  actual  temperature  at  the  lower 

entrance  and  convert  to  virtual  temperature  using 
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the  specific  humidity  for  atmospheric  conditions. 

6)  Calculate  the  initial  Mach  number. 

7)  Calculate  the  initial  pressure  and  density. 

8)  Calculate  the  mass  rate  of  flow  and  the  relaxation 
length. 

9)  Begin  calculating  parameters  of  subsequent  cross- 
sections  by  calculating  for  each  cross-section  the 
specific  humidity,  the  saturated  specific  humidity, 
and  the  temperature  according  to  the  appropriate 
formulas. 

10)  Calculate  the  virtual  temperature  and  iterate  for 
the  mach  number  of  each  cross-section. 

11)  Calculate  pressure  at  each  cross-section. 

12)  Continue  until  end  of  cave  is  reached  -  Is  the  final 

pressure  equal  to  atmospheric  pressure?  -  If  no, 
estimate  new  initial  velocity  from  the  temperature 
and  velocity  calculated  and  the  actual  atmospheric 
pressure  -  Return  to  step  4  -  If  yes,  calculate  all 
remaining  parameters  from  known  Mach  number 
humidity,  and  temperature  variations. 

Case  studies  of  this  and  other  models  are  contained  in 


section  VI. 
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V  -  INCOMPRESSIBLE  MODEL 


The  compressibility  of  air  severely  complicates  the 
determination  of  speleomicroclimatic  parameters  as  is 
witnessed  by  the  foregoing  model.  It  would  be  of  great  value 
if  air  could  be  treated  as  imcompressible  for  the  purpose  of 
calculating  these  parameters.  This  is  not  an  unwarranted 
assumption  for  level  passages  that  are  short  relative  to 
their  diameter  (the  effect  of  friction  on  the  flow  is  not 
extremely  strong)  if  the  mach  number  is  small.  Shapiro 
(1953,  p.48)  states  that  "...the  error  produced  by 
neglecting  compressibility  in  the  computation  of  pressure 
variations  is  of  the  order  of  one-fourth  the  square  of  the 
ratio  of  the  stream  velocity  to  the  sound  velocity  [Mach 
number]."  If  the  given  requirements  were  met,  the  assumption 
of  incompressibility  should  produce  an  error  in  the 
calculation  of  pressure  of  roughly  only  one  per  cent  at  a 
Mach  number  of  0.2.  An  error  of  this  magnitude  is  more  than 
acceptable  for  the  results  desired  here,  but  unfortunately 
the  assumptions  are  not  met.  With  a  level  passage,  no 
consistent  chimney  effect  can  exist,  and  most  caves  which 
exhibit  this  effect  are  generally  quite  long  with  respect  to 
the  size  of  the  passage.  The  range  of  Mach  numbers  that  are 
dealt  with  here,  however,  are  so  low  (0.0005  -  0.05)  that 
these  may  not  be  sufficient  objections,  so  an  incompressible 
model  is  worth  investigating. 

It  is  obvious  that  the  density  of  air  must  vary  with 
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elevation  in  a  non- horizontal  passage,  so  the  air  cannot  be 
treated  as  totally  incompressible.  The  assumption  being  made 
is  that  dynamic  influences  on  density  are  not  significant. 
It  is  necessary,  then,  to  know  how  density,  or  equivalently 
pressure,  varies  with  elevation  under  static  conditions. 
Given  no  flew  of  air,  there  can  be  no  temperature 
differential  between  the  cave  air  and  the  surrounding  rock, 
so  the  internal  atmosphere  is  at  the  internal  wall 
temperature,  TWALL.  Since  TW ALL  has  been  assumed  constant 
throughout  the  cave,  the  cave  atmosphere  (under  the 
assumption  cf  no  flow  only)  is  isothermal,  and  the  pressure 
variation  with  height  may  be  calculated  from  the  hydrostatic 
equation 

dP  =  :  -  DGdz  .  (5.1) 

Replacing  D  by  P/(Rg)T  from  equation  (4.1)  and  rearranging, 
we  get 

dP  =  : -Gdz_  , 

P  "(Eg)  T  (5.2) 

which  may  be  integrated  to  yield 

PJ2]_  =  exp  (-G  (z  (2)  -  z  (1)  )/(  (Rg)  TWALL) ) 

P(1)  (5-3) 

Were  the  temperature  variation  independent  of  velocity,  the 
assumption  cf  incompressibility  would  imply  that  this 
equation  holds  for  dynamic  as  well  as  static  conditions.  The 

variation  is  strongly  dependent  upon  the  air 
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velocity/  however/  and  hence  density  is  indirectly  dependent 
upon  this  velocity.  The  assumption  of  incompressibility 
must/  in  this  situation,  mean  that  the  actual  and  stagnation 
densities  are  the  same,  not  that  the  density  is  invariant. 
The  situation  may  be  more  clearly  understood  by  picturing 
the  cave  atmosphere  as  being  composed  of  horizontal  layers 
of  infinitesimal  thickness.  Air  density  may  vary  from  layer 
to  layer,  but  is  constant  within  any  given  layer  (Thompson/ 
1972) .  The  effect  of  a  change  in  flow  conditions,  with  a 
resultant  change  in  the  temperature  distribution,  serves  not 
to  change  the  actual  densities,  but  merely  to  shift  the 
layers  upward  or  downward.  It  is  this  upward  and  downward 
motion  (of  the  density  layers)  that  allows  the  cave 
atmosphere  to  adjust  to  external  conditions.  The  solution 
for  the  density  variation,  then,  is  very  similar  to  that 
under  static  conditions,  but  the  temperature  variation  is 
now  that  of  equations  (4.58)  ,  (4.61)  ,  (4.62)  ,  and  (4.63)  • 

The  solution  from  equations  (4.30)  and  (5.2)  requires 
numeric  integration  (trapezoidal  rule) ,  since  T  is  no  longer 
constant.  This  yields 

r  r  t  r  n 

D  (2)  =  D  (1)  T  (2)  exp|  -GsinO  I  _J. _ +  _J _  II  2Lp2j_-_x j  1j_  |  j  , 

T  ( 1)  |  (Hg)  I  1(2)  T  (1)  ||  2  || 

l  t  J  L  JJ 

(5.4) 

inhere,  as  before,  z  =  sin6x,  so  dz  —  sinOdx.  Given,  then, 
the  temperature  variation  along  the  passage,  all  that 
remains  now  is  to  interpret  the  boundary  conditions  in  terms 
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o±  the  flow  parameters.  In  the  compressible  model,  the 
atmospheric  conditions  were  taken  as  the  initial  stagnation 
conditions,  so  let  us  look  at  equations  (4.25)  and  (4.26), 
which  are 


To 


+  JkrJIM  I  , 
2  l 


and 


(5.5) 


r  nk/(k-1) 

Po  =  P|  1  ♦  Jk^JlMZ  | 

I  2  | 

1  J  (5.6) 

The  highest  Mach  number  that  we  would  be  concerned  with  in 
cave  studies  is  around  0.05,  which  corresponds,  at  a 

temperature  of  10°C,  to  an  air  velocity  of  approximately  17 
meters/sec.  Since  T  and  To  are  measured  in  °K,  we  may 
generously  assume  that  the  range  of  temperatures  cf  interest 
to  us  is  230-310°K,  so  for  To  varying  within  this  range,  the 
maximum  error  in  neglecting  M  is  cnly  0.11  to  0.15°K,  a 
minimal  difference.  This  result  was  anticipated  in  the 
development  of  equations  (5.3)  and  (5.4),  since  it  was 
assumed  without  explanation  that  the  temperature  did  not 
vary  due  to  compressibility  of  the  air.  For  P,  however,  the 
power  term  k/(k-1)  is  equal  to  3.5,  and  the  dependence  upon 
M  cannot  be  ignored.  Instead,  let  us  binomially  expand  the 
right  side  of  equation  (5.6) .  This  gives 
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Po 


r  ! 

=  P|  1  +  kM£  +  kfl*  +  k (2-k) M6  + 

I  2  8  ”  ~48  “ 

L  j 

=  F  +  PkMg  ♦  Pkj4  +  Pk(2-k) ♦  ... 

2  8  48 


=  P  +  PkM£| 

2  I 


♦  M  +  (2-k)M*  ♦  ...  ]  . 

2  24  j 
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But  from  eguation  (4.3), 


(5.7) 


PkM£  =  _PkV£ _ 

2  2k  (Eg)  T 

=  _PV£ _ , 

2  (Rg)  T 

and  since 

D  =  P _ 

(Bg)T 


from  equation  (4.1),  it  follows  that  equation  (5.7)  becomes 


Po 


=  P  ♦ 


DV£| 
2  I 


1  + 


M 

4 


(2-k)  M4 

24 


(5.8) 


Neglecting  the  Mach  number  in  this  equation  induces  no  great 
error,  hence  we  arrive  at 


Po  =  P  +  DV£  , 

2  (5.9) 

which  is  familiar  to  all  those  acquainted  with  fluid 
dynamics  as  the  expression  for  the  stagnation  pressure  for 
incompressible  flow.  P  is  often  called  the  static  pressure 
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and  DV2/2  the  velocity  pressure.  The  stagnation  pressure,  in 
this  context,  is  often  called  the  total  pressure  for  obvious 
reasons. 

The  solution  is  now  available  from  the  Eernoulli  - 
Euler  equation 

0  =  dP  +  DVdV  +  DGsin-edx  ♦  fDV2dx 

R  (5.10) 

by  integration  over  the  length  of  the  cave.  This  yields 
(using  Simpson^  one-third  rule), 

0  =  P(N  +  1)  -  P  (1)  +  D(1)V(1)  (V(2)  -  V(1)) 

+  GsinO  (D  (1)  ♦  4D  (2)  *  ...  +  D(N  +  1))(x(2)  -  X(1))/3 
+  fV  (1)  ((1/0(1))  +  (4/D  (2)  )  *  ...  ♦  (1/D  (N  + 1)  ) 

•  (X  (2)  -  X  (1 )  )  /  (3R)  (5.11) 

since  DV  is  constant  when  w  and  R  are  constant.  Rearranging 
and  noting  that 

P(1)  =  PEXTL  -  0.5D  (1)  V  (I)2  , 

we  get 

P  (N  +  1) 

=  PEXTL  -  0  •  5D  ( 1)  V(1)  2  -  D  (1)  V  (1)  2  (  (D  (1)/D  (2)  )  -  1) 

-  GsinO  (D  (1)  ♦  4D  (2)  ♦  ...  *  D(M*1))(x(2)  -  X(1))/3 

-  fV  (1)  (  (1/D  (1)  )  +  (4/D  (2 )  )  +  ...  ♦  (  1/D  (N  +  1)  ) 

•  (x  (2)  -  x(1))/(3R)  (5.12) 


From  this. 


the  initial  velocity  V  (1 )  may  be  derived  by 


. 


. 

. 


1 

. 


... 


■ 

■ 


.  • 

•  '  1 

I 


56 


iteration,  comparing  P(N+1)  to  PEXTU  until  they  are 
f icient ly  close  (so  that  each  calculated  velocity  is 
within  0.0001  m/sec  of  the  previously  calculated  velocity), 
tfith  the  temperature  variation  obtained  from  eguations 
(4.58),  (4.61),  (4.62),  and  (4.63),  the  solution  is  now 
readily  available  from  equations  (5.4)  and  (5.9),  and  the 
equations  fcr  temperature  distribution,  (4.58),  (4.61), 
(4.62)  and  (4.63) .  The  procedure  is  as  follows 

1)  Specify  the  external  atmospheric  conditions. 

2)  Specify  the  cave  parameters. 

3)  Estimate  the  initial  velocity  and  calculate  the 
temperature  and  hence  pressure  variation  step  by 
step  along  the  passage. 

4)  From  the  external  pressure  and  calculated  air 
temperature  at  the  upper  entrance,  calculate  the 
pressure  at  the  lower  entrance. 

5)  Compare  the  calculated  pressure  at  the  lower 
entrance  with  the  stagnation  pressure  (atmospheric) 
there,  and  estimate  the  new  initial  velocity. 

6)  Continue  the  iteration  on  the  initial  velocity  until 
the  desired  accuracy  is  reached. 

7)  Calculate  all  remaining  parameters. 

The  above  procedure  has  been  simplified  by  the  ignoring  of 
moisture,  and  the  reader  who  wishes  a  more  complete 
procedural  description  should  compare  the  above  with  that 
for  the  compressible  model  on  pages  48  and  49. 
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Unconstrained  Temperature  Profile 

A  further  simplification  which  can  be  made  to  the 
incompressible  model  without  severe  loss  of  accuracy  is 
suggested  by  Wigley  and  Brown  (1971).  Removal  of  the 
constraint  that  the  air  may  not  become  saturated  produces 
only  minor  error  and  allows  the  use  of  equations  (4.58)  and 
(4.62)  for  the  temperature  variation  along  the  entire 
passage,  thus  avoiding  the  iterative  solution  of  equation 
(4.61).  "[One  is]  justified  in  using  the  non-constrained 
solution  for  the  temperature  distribution  in  the  majority  of 
quantitative  applications"  (Wigley  and  Brown,  1971,  p.312). 
The  procedure  here  is  identical  to  that  in  the  previous 
section.  Equation  (5.12)  again  used  to  iterate  for  the 
initial  velocity,  with  the  temperature  variation  coming  from 
equations  (4.58)  and  (4.62),  depending  on  whether  or  not  the 
effects  of  moisture  are  considered.  A  comparison  of  this  and 
the  previous  models  is  found  in  the  following  section. 
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VI  -  CASE  STUDIES 


In  order  to  compare  the  three  previously  described 
models,  the  following  hypothetical  situation  is  constructed. 
Assume  that  the  ideal  dynamic  cave  described  in  section  IV 
is  1000  metres  long  and  rises  at  a  constant  angle  of  10°. 
The  internal  wall  temperature  is  10°C,  the  wall  friction 
factor  is  0.030,  and  the  radius  of  the  passage  is  1.5 
metres.  The  walls  of  the  passage  will  be  considered  wet 
throughout  its  length,  so  E  =  1.0  .  In  order  to  determine 
the  internal  microclimatic  conditions,  all  that  is  required 
are  the  external  atmospheric  conditions.  Three  cases  were 
devised  which,  for  simplicity  in  the  comparison,  vary  only 
in  external  temperature.  These  three  cases  are: 

Atmospheric  pressure  at  upper  entrance  =  925  mb 

Atmospheric  pressure  at  lower  entrance  =  950  mb 

Relative  Humidity  at  lower  entrance  =  0.50 

Temperature  at  lower  entrance  (a)  =  -20  °C 

(b)  =  -10  °C 

(c)  =  0  °C 

The  internal  parameters  were  calculated  at  25  metre 
intervals  and,  for  simplicity,  the  variation  of  the  specific 
humidity  (at  saturation)  with  pressure  was  ignored.  The 
results  are  tabulated  in  Appendix  A,  while  the  computer 
programs  used  are  found  in  Appendix  B. 
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The  Compressible  Model 

As  this  study  is  entirely  theoretical  in  nature,  the 
norm  by  which  models  are  accepted  or  rejected  must  of 
necessity  be  that  one  which  makes  the  fewest  assumptions 
about  the  nature  of  the  flow,  namely  the  compressible  model. 
The  three  cases  studied  all  show  very  little  variation 
between  density,  pressure,  and  temperature  and  their  related 
stagnation  properties.  The  difference  is  on  the  order  of 
0.05  per  cent  of  the  stagnaticn  properties  and  hence  may 
generally  be  ignored.  The  reason  for  this  minute  difference 
is  evident  from  the  Mach  numbers  calculated.  In  no  instance 
is  the  Mach  number  greater  than  0.015,  and  the  variation 
along  the  passage  is  minimal.  Thus  the  relationship  between 
velocity  and  the  state  parameters  may  be  ignored  (with  the 
obvious  exception  of  temperature,  which  depends  on  the  rate 
of  heat  transferral  to  or  from  the  air,  a  velocity  dependent 
operation) . 

The  relative  humidity  bears  special  note,  for  it  almost 
immediately  reaches  100  per  cent.  This  is  to  be  expected, 
however,  for  under  winter  conditions  such  as  those  here, 
there  is  very  little  water  vapor  in  the  air  at  saturation. 
Upon  entering  the  cave,  the  air  very  quickly  gains  moisture 
from  the  cave  walls  and  saturation  is  achieved  almost  at 
once.  The  opposite  effect  occurs  in  summer  for  the  upper 
entrance,  where  the  air  is  capable  of  holding  considerably 
more  water  vapor  than  under  winter  conditions.  Upon  entering 
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the  cave,  the  air  takes  on  moisture  from  the  vails  for  a 
much  longer  period  of  time  (than  under  winter  conditions 
with  the  same  relative  humidity  at  the  entrance  drawing  air) 
before  saturation  is  finally  (if  ever)  reached  (if  the 
external  air  is  saturated ,  however,  then  the  entire  cave 
will  also  be  saturated,  regardless  of  season) . 

By  far  the  most  curious  result  obtained  from  the 
compressible  model,  however,  is  the  air  velocity.  While  the 
overall  magnitude  of  the  velocity  in  each  case  is  reasonable 
(in  light  of  the  observed  magnitude  of  most  cave  winds) 
considering  the  simplif ications  imposed  (no  constrictions, 
no  bends,  etc.),  the  relative  magnitude  between  cases 
appears  contradictory.  As  the  external  temperature 
approaches  the  internal  wall  temperature,  the  internal- 
external  pressure  differential  decreases,  so  the  air 
velocity  should  also  decrease.  Exactly  the  opposite, 
however,  is  found  in  the  three  cases  studied;  as  the 
external  temperature  increases  from  — 20°C  to  -10°C,  the 
initial  velocity  increases  from  4.25  m/sec  to  4.42  m/sec  to 
4. 57m/sec.  The  contradiction  is  only  apparent,  however,  and 
is  due  entirely  to  the  cases  chosen.  Because  the  external 
atmospheric  pressure  at  both  entrances  was  specified  as 
constant,  an  increase  in  temperature  at  the  lower  entrance 
must  bring  about  a  decrease  in  air  density  there.  Thus,  to 
maintain  the  designated  pressure  differential  between  the 
entrances,  the  less  dense  air  must  move  faster  to  have  the 


equivalent  effect.  The  reader  should  thus  be  aware  that  the 
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cases  under  study  are  not  meant  to  be  related  in  any  manner 
through  time*  but  rather  are  distinct  situations  in  which 
certain  variables  are  held  constant  so  that  the  others  may 
be  more  fully  understood,  in  an  actual  situation*  one  would 
expect  seme  variation  in  pressure  to  accompany  such  a 
variation  in  temperature. 

The  Incompr e s si ble  Models 

Comparison  of  the  results  generated  by  the 
incompressible  models  with  those  produced  by  the 
compressible  model  shows  excellent  agreement  with  regard  to 
pressure  variations.  The  difference  between  all  models  in 
this  respect*  is  again  on  the  order  of  0.02  per  cent*  but 
this  could  easily  be  anticipated.  The  pressure  differential 
between  the  entrances  is  only  25  millibars  and  furthermore* 
the  entrance  pressures  are  fixed  limits  between  which  each 
model  must  operate.  Considering  that  the  printed  results  are 
accurate  only  to  0.1  millibars*  one  would  readily  expect  the 
pressure  variation  along  the  passage  to  differ  little  from  a 
mean  of  2.5  millibars  per  hundred  metres*  and  such  is  in 
fact  the  case.  The  variation  over  all  cases  and  all  models 
is  2.2  to  2.6  with  the  variation  consistently  greater  at  the 
lower  entrance  and  less  at  the  upper.  The  predictive  ability 
of  the  incompressible  models*  then*  cannot  be  judged  on 
their  ability  to  predict  pressure  variations. 


When 


temperature  and  density  are  considered* 
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differences  between  the  models  begin  to  appear.  The 
constrained  incompressible  model  still  agrees  well  with  the 
compressible ,  the  variation  of  the  former  over  the  latter 
being  on  the  order  of  0.2  per  cent.  The  unconstrained 
incompressible  model,  however,  varies  nearly  3.0  per  cent 
from  the  compressible  model  in  both  temperature  and  density. 
While  errors  of  this  magnitude  might  not  appear  to  be 
significant  in  light  of  the  numerous  simplifications  used  to 
generate  them,  it  should  be  noted  that  for  the  analysis  of 
error,  temperatures  must  be  expressed  in  degrees  Kelvin  and 
a  3.0  per  cent  error  at  the  temperatures  under  consideration 
is  approximately  8.0°C.  This  difference  is  most  readily  seen 
in  case  1,  where  the  temperature  at  150  metres  is  calculated 
at  -8.46°C  in  the  unconstrained  incompressible  model  and  at 
-1.44°C  in  the  compressible  model.  While  Wigley  and  Brown 
(1971)  have  shown  that,  at  the  same  velocity,  and  with  a 
relatively  small  temperature  differential  between  the 
atmosphere  and  the  cave  wall,  the  constrained  and 
unconstrained  temperature  profiles  are  sufficiently  similar, 
here  the  temperature  differential  is  extreme  (30°C)  and 
markedly  erroneous  estimations  of  temperature  are  produced. 

In  general,  the  constrained  incompressible  model  agrees 
well  with  the  compressible  model.  Temperature,  pressure,  and 
density  vary  less  than  0.5  per  cent  from  the  compressible 
values,  while  velocity  varies  on  the  order  of  3.2  per  cent. 
In  nearly  all  applications,  then,  the  constrained 
incompressible  model  presents  an  excellent  approximation  to 
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the  compressible. 

A  word  of  caution  concerning  velocity  estimations  from 
the  above  models  is  in  order.  Wigley  and  Brown  (1971)  have 
pointed  out  that  the  relaxation  length  for  temperature  is 
very  insensitive  to  changes  in  velocity  (the  dependence  is 
only  to  the  0.2  power) .  The  pressure  variation  is  also 
weakly  dependent  upon  velocity,  for  it  was  noted  during  the 
iterations  for  initial  velocity  that  a  reduction  of  only  two 
millibars  in  the  calculated  pressure  at  the  upper  entrance 
required  an  increase  in  velocity  of  nearly  1.0  metres  per 
second.  Very  accurate  measurements  or  estimates  of  the 
pressure  differential  between  the  entrances  would  thus  be 
necessary  for  the  reasonable  estimation  of  velocities  within 
the  cave. 

The  results  of  the  three  case  studies  are  displayed 
graphically  in  figures  11  to  20,  and  the  actual  computer 
print-out  and  the  programs  used  appear  in  Appendices  A  and  B 
respectively . 
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VII  -  CONCLUSION 


Context  and  Contribution  of  Present  Study 

It  has  been  shown  that  the  constrained  incompressible 
model  is  an  excellent  approximation  to  the  compressible 
model  for  the  flow  of  air  in  caves.  At  no  point  in  the  case 
studies  did  the  two  differ  beyond  the  bounds  of  the  error 
inherent  in  the  numerical  solutions,  and  there  is  no  reason 
to  expect  any  greater  difference  if  these  solutions  were 
refined.  Despite  the  long  distances  over  which  friction 
works,  and  despite  the  substantial  changes  in  elevation,  it 
is  clear  that  the  mach  number  is  sufficiently  low  to  allow 
use  of  a  modified  incompressible  formulation.  The  Bernoulli- 
Euler  equation  in  its  compressible  form,  the  hydrostatic 
equation,  and  the  temperature  variation  derived  by  Wigley , 
and  Brown  (1971)  give  a  more  than  adequate  determination  of 
the  microclimatic  variables  in  our  ideal  cave. 

Any  extension  of  this  model  to  actual  caves  must  be 
done  with  caution,  however.  The  reader  is  reminded  that  the 
present  study  is  an  entirely  theoretical  development  in 
which  numerous  simplifying  assumptions  have  been  made.  In 
relatively  short,  unconstricted ,  non*sinuous  caves  with 
little  vertical  extent,  the  present  theory  may  be  expected 
to  yield  reasonable  results.  If  these  conditions  are  not 
met,  the  agreement  of  the  model  with  reality  may  be  poor. 
The  effects  of  sinuosity  may  be  accommodated  by  the  friction 
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factor,  for  the  effect  of  bends  in  the  passage  is  to  resist 
the  flow  of  air.  Changes  in  cross-sectional  area,  however, 
would  be  somewhat  more  bothersome.  Their  inclusion  into  the 
compressible  model's  equation  for  mach  number  variation 
(4.32)  is  straightforward.  Wigley  and  Brown  (1971),  however, 
have  pointed  out  that  their  inclusion  into  the  equations  for 
temperature  variation  makes  the  relaxation  length  variable 
and  may  negate  the  empirical  relationship  between  the 
Nusselt,  Prandtl,  and  Reynolds  numbers.  Therefore,  given  the 
present  state  of  heat  and  mass  transfer  theory,  an 
assumption  of  constant  radius  (i.e.,  an  average  value)  may 
be  necessary. 

In  actual  caves,  changes  in  elevation  and,  more 
important,  changes  in  depth  beneath  the  surface  produce 
changes  in  the  rock  temperature,  making  TWALL  non-constant. 
The  requisite  differential  equations  may  be  integrated 
numerically  assuming  a  non— constant  wall  temperature,  but 
again  the  empirical  relationship  between  the  Nusselt, 
Prandtl,  and  Reynolds  numbers  may  not  hold.  From  the  study 
of  Kays  (1966),  it  would  appear  that  this  relationship  may 
still  be  a  reasonable  approximation,  but  a  more  thorough 
knowledge  of  heat  and  mass  transfer  theory  than  the  author 
possesses  would  be  necessary  to  pursue  this  point. 

In  sum,  then,  this  thesis  analyses  the  problem  and  the 
constraints  involved,  rather  than  presenting  a  solution  as 
such  to  that  problem.  The  modified  incompressible  model 
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presented  is  only  a  basic  framework.  Refinements  to  the 
theory  are  needed  but,  more  important,  substantial  fieldwork 
is  called  for  to  adapt  this  model  to  reality.  Only  through 
well* planned  empirical  studies  may  the  theory  more  closely 
approximate  actuality. 

future  Research 

Concurrent  with  the  wide  acceptance  of  static  theories 
for  cave  air  flow  is  a  lack  of  carefully  executed  fieldwork. 
The  author  knows  of  no  instance  where  the  air  temperature, 
pressure,  and  velocity  within  a  cave  have  been  measured 
continuously  over  a  long  period  of  time  at  more  than  one 
location.  The  complete  monitoring  of  a  cave  at  numerous 
locations  along  its  length  needs  to  be  undertaken,  but  there 
are  few  caves  in  Canada  that  are  sufficient  for  such 
studies.  Nakimu  Caves  in  Glacier  National  Park,  British 
Columbia,  may  be  adequate,  but  the  presence  of  a  substantial 
stream  in  the  cave  complicates  the  temperature  profile  and 
alters  the  air  velocity  through  entrainment.  Grotte  Valerie, 
Northwest  Territories,  has  also  been  suggested,  but  it  is 
over  70  miles  from  the  nearest  highway  and  thus  presents 
considerable  problems  in  logistics.  Certain  caves  on 
Vancouver  Island,  British  Columbia,  may  also  suffice,  though 
they  tend  to  be  too  short  or  too  constricted.  Nevertheless, 
the  behavior  of  these  caves  under  changing  external 
conditions  needs  to  be  monitored  carefully  to  shed  light  on 
the  empirical  relationships  between  the  numerous  parameters 
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involved . 

The  extension  of  the  steady-flow  case  to  unsteady 
situations  needs  also  to  be  undertaken.  Wigley  and  Brown 
(1971)  have  attempted  to  predict  temperature  profiles  under 
changing  external  conditions,  but  they  assume  a  constant 
velocity,  which  can  at  best  be  considered  only  an 
approximation.  The  solution  of  the  unsteady  case  is  of  great 
importance  for  the  use  of  airflow  in  caves  in  the  remote 
sensing  of  subsurface  karst  characteristics,  as  well  as  in 
the  determination  of  variability  of  the  subterranean 
atmosphere.  It  is  this  combination,  then,  theory  and 
empiricism,  that  will  ultimately  lead  to  a  better 
understanding  of  subterranean  microclimatic  conditions. 
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Table  Is 

SYMB01ISATI0N  AND 

UNITS 

S  ymbol 

Units 

Meaninq 

A 

m2 

cross-sectional  area  of 

passage 

A£  wall  ] 

m2 

area  of  passage  wall 

B 

m 

hydraulic  diameter 

c£p] 

m2/sec2-°K 

specific  heat  of  dry  air  at 
constant  pressure 

D 

kg/m3 

mass  density 

Do 

kg/m3 

stagnation  density 

dh[  pr] 

m2/sec2 

change  in  enthalpy  for  change 
in  state 

E 

- 

fraction  of  passage  wall  that 
is  moist 

£ 

- 

Fanning  friction  factor 

f £  smooth  ] 

- 

friction  factor  for  smooth 

passage 

G 

m/sec2 

acceleration  due  to  gravity 

H 

m2/sec2 

latent  heat  of  evaporation  (or 
sublimation) 

h 

m2/sec2 

specific  enthalpy 

ho 

m2/sec2 

stagnation  enthalpy 

k 

• 

ratio  of  specific  heat 
constant  pressure  to  specific 
heat  constant  volume 

L 

°K/m 

dry  adiabatic  lapse  rate 

m 

- 

mach  number  -  fraction  of 

speed  of  sound 

Nu 

- 

Nusselt  number 

p 

kg/m -sec2 

pressure 

Po 

kg/m-sec2 

stagnation  pressure 
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PEXTL 

kg/m-sec2 

PEXTU 

kg/m-sec2 

Pr 

- 

Q 

m2/sec2 

q 

- 

QEXT 

- 

QSAT[T] 

- 

QSPEC 

- 

QW  ALL 

- 

R 

m 

Re 

- 

Rg 

m2/sec2-°K 

fiv 

m2/sec2- °K 

T 

°K 

t 

sec 

t[ wall  ] 

kg/m-sec2 

To 

°K 

Tp 

°K 

TEXTL 

°K 

TSAT 

°K 

Tv 

°K 

THALL 

°K 

V 

m/sec 

w 

kg/sec 

atmospheric  pressure  at  lower 
entrance 

atmospheric  pressure  at  upper 
entrance 

Prandtl  number 

heat  per  unit  mass  of  air 

specific  humidity 

specific  humidity  at  lower 
entrance 

specific  humidity  at 

saturation  for  temperature  T 

mean  specific  humidity  of  cave 
air 

specific  humidity  at 

saturation  for  temperature 
TWALL 

radius  of  passage 

Reynolds  number 

gas  constant  for  dry  air 

gas  constant  for  water  vapor 

temperature 

time 

wall  shear  stress 
stagnation  temperature 
potential  temperature 
temperature  at  lower  entrance 
temperature  at  XSAT 
virtual  temperature 
wall  temperature  in  cave 
air  velocity 
mass  rate  of  flow 
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hrs. 
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w 


wavelength  of  temperature 
variation 


X 


x 


kg-m/sec*  component  of  gravitational 

force  working  against  flow 

m  distance  into  cave 


xo 


xsat 


m  relaxation  length  -  distance 

at  which  temperature 

differential  has  decayed  to 
1/e  its  original  value 

m  distance  into  cave  at  which 

saturation  is  achieved 


Z[x] 


m 


elevation  of  point  in  cave  as 
a  function  of  x 


Z  m 

0 


elevation  above  lower 

entrance;  depth  below  surface 

angle  of  inclination  of 
passage 

denotes  value  of  variable 
corresponding  to  cross-section 
at  distance  x (i)  into  cave 


(i) 
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FIGURE  1:  Idealised  Static  Cave 

P  [ext]  >  P  [int] 


P  [ext]  <  P  ['nt] 


V 


J  . 
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FIGURE  2:  Dynamic  Cave  in  Equilibrium  with  Atmosphere 

(after  Cigna,  1963) 
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FIGURE  3:  Air  Columns  of  Different  Weights 

-  Summer  Conditions  in  a  Dynamic  Cave 
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FIGURE  4:  Mine  with  Artificial  Ventilation -Winter  Conditions 

(Enhanced  Natural  Ventilation) 
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FIGURE  5:  IDEAL  PRESSURE- SPECIFIC  VOLUME  GRAPH 
(THE  JOULE  OR  BRAYTON  CYCLE  (AFTER 
HALL,  1967) 
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FIGURE  6:  Actual  Pressure  -  Specific  Volume  Graph 
(after  Hall,  1967) 
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FIGURE  7 
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:  Pressure  -  Specific  Volume  Graph  for  Step- 

by -  Step  Analysis  (Gfter  Hall,  1967) 
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FIGURE  8:  Ideal  Dynamic  Cave 

(see  Table  1  for  explanation  of  symbols) 
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FIGURE  9:  The  Control  Volume  (after  Shapiro,  1953) 


' 


80 


FIGURE  10:  Approximation  to  Geothermal  Gradient 
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FIGURE  II  =  Temperature  -  Case 
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FIGURE  12  :  Temperature  -  Case 
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FIGURE  13  •*  Temperature  -  Case 
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FIGURE  14  :  Velocity  -  Case 


84 


r— — i - r 

O  CD 

to  <3- 


1 - 1 - 1 - 1 —  I 

tO  ^  CM 

'T  ^ 


"1 

O 


(09S/UJ)  A  1 1  3  0 1 3A 


O  compressible 
A  constrained 
□  unconstrained 


85 


I 


OJ 


CD 

to 

O 

O 


>» 

4— 

o 

o 

Q) 

> 


in 


LlJ 

cr 

3 

O 

Lu 


o 

o 

OJ 


o 

o 


LU 

O 

2 

< 

h- 

co 

o 


O  <3  □ 


i - 1 - r 

O  00 

lO 


1 - 1 - 1 - 1 - 1 - 1  - 1 

<X)  ^  o 


(  033  /  LU  ) 


A.1IOOT3A 


compressible 

constroined 

unconstrained 


FIGURE  16=  VELOCITY  -  Case 
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FIGURE  17  :  Density  -  Case 
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FIGURE  18  :  Density  -  Case 
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FIGURE  19**  Density  -  Case 
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FIGURE  20  :  Pressure 
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APPENDIX  A:  CASE  STUDIES 


Case  1;  Compressible  mode 1 


TEMPERATURE  AT  LOWER  ENTRANCE 
PRESSURE  AT  LOWER  ENTRANCE 
PRESSURE  AT  UPPER  ENTRANCE 
RELATIVE  HUMIDITY  AT  LOWER  ENTRANCE 
NUMBER  OF  DIVISIONS  OF  PASSAGE 


-20.000  DEG  C. 
950.000  MB 
925.000  MB 
0.500 
40 


LENGTH  OF  PASSAGE 
RADIUS  OF  PASSAGE 
ANGLE  OF  PASSAGE 
WALL  TEMPERATURE 
FRICTION  FACTOR 

FRACTION  OF  WALL  THAT  IS  MOIST 


1000.000  M 
1.500  M 
10.000  DEG 
10.000  DEG  C. 
0.030 
1.000 


DISTANCE  TEMPERATURE  PRESSURE  DENSITY 

(METRES)  (DEG  CELSIUS)  (MILLIBARS)  (KG/CUBIC  M) 


0.0 

-19.99 

25.0 

-15.52 

50.0 

-11.7 2 

75.0 

-8.49 

100.0 

-5.75 

125.0 

-3.42 

150.0 

-1.44 

175.0 

0.24 

200.0 

1.67 

225.0 

2.89 

250.0 

3.92 

275.0 

4.79 

300.0 

5.54 

325.0 

6.17 

350.0 

6.71 

375.0 

7.  16 

400.0 

7.55 

425.0 

7.88 

450.0 

8.16 

4  75.0 

8.40 

500.0 

8.60 

525.  0 

8.77 

550.0 

8.92 

575.0 

9.04 

600.0 

9.15 

625.  0 

9.24 

650.0 

9.31 

675.0 

9.38 

700.0 

9.43 

725.0 

9.48 

750.0 

9.52 

775.0 

9.55 

800.0 

9.58 

825.0 

9.60 

950.0 

1.3073 

949.3 

1.2836 

948.7 

1.2642 

948.1 

1.2479 

947.5 

1.2342 

946.9 

1.2227 

946.2 

1.2130 

945.6 

1.2047 

945.0 

1.1976 

944.3 

1.1915 

943.7 

1.1862 

943.1 

1.1817 

942.5 

1.1777 

941.8 

1.1742 

941.2 

1.1712 

940.6 

1.1685 

939.9 

1.1661 

939.3 

1.1639 

938.7 

1.1620 

938.1 

1.1602 

937.4 

1.1586 

936.8 

1.1571 

936.2 

1.1557 

935.6 

1.  1545 

934.9 

1.1532 

934.3 

1.1521 

933.7 

1.1510 

933.1 

1.1500 

932.4 

1.1490 

931  .8 

1 .1480 

931.2 

1.1471 

930.6 

1. 1462 

930.0 

1.1453 

929.3 

1 .1445 
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850.  0 

9.62 

875.0 

9.64 

900.0 

9.66 

925.0 

9.67 

950.0 

9.68 

975.0 

9.68 

1000.0 

9.70 

STAGNATION 

DISTANCE 

TEMPERATURE 

(METRES) 

(DEG  CELSIUS) 

0.0 

-19.98 

25.0 

-15.51 

50.0 

-11.71 

75.0 

-8.48 

100.0 

-5.74 

125.0 

-3.41 

150.0 

-1.43 

175.0 

0.25 

200.0 

1.68 

225.0 

2.90 

250.0 

3.93 

275.0 

4.80 

300.0 

5.55 

325.0 

6.18 

350.0 

6.72 

375.0 

7.17 

400.0 

7.56 

425.0 

7.89 

450.0 

8.  17 

475.0 

8.41 

500.0 

8.6  1 

525.0 

8.78 

550.0 

8.93 

575.0 

9.05 

600.0 

9.16 

625.0 

9.25 

650.0 

9.32 

675.0 

9.39 

700.0 

9.44 

725.0 

9.49 

750.0 

9.53 

775.0 

9.56 

800.0 

9.59 

825.0 

9.61 

850.0 

9.64 

875.0 

9.65 

900.0 

9.67 

925.0 

9.68 

950.0 

9.69 

975.0 

9.69 

1000.0 

9.71 

928.7 

1.1436 

928.1 

1.1428 

927.5 

1.1420 

926.9 

1.1411 

926.2 

1.1403 

925.6 

1.1396 

925.0 

1.1387 

STAGNATION 

STAGNATION 

PRESSURE 

DENSITY 

(MILLIBARS) 

(KG/CUBIC  M) 

950.1 

1.3074 

949.5 

1.2838 

948.9 

1.2643 

948.2 

1.2480 

947.6 

1.2344 

947.0 

1.2229 

946.4 

1.2131 

945.7 

1.2048 

945.1 

1.1977 

944.5 

1.1916 

943.8 

1.1863 

943.2 

1.1818 

942.6 

1.1778 

942.0 

1.1744 

941  .3 

1.1713 

940.7 

1.1686 

940.1 

1.1662 

939.5 

1.1640 

938.8 

1.1621 

938.2 

1.1603 

937.6 

1.1587 

936.9 

1.1572 

936.3 

1.1559 

935.7 

1.1546 

935.1 

1.1534 

934.5 

1.1522 

933.8 

1.1511 

933.2 

1.1501 

932.6 

1.1491 

932.0 

1.1482 

931  .3 

1.1472 

930.7 

1.1463 

930.1 

1.1454 

929.5 

1.1446 

928.9 

1. 1437 

928.2 

1.1429 

927.6 

1.1421 

927.0 

1.1413 

926  .4 

1.1404 

925.7 

1.1397 

925.1 

1.1389 
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DISTANCE 

VELOCITY 

(METRES) 

(M/SEC) 

0.0 

4.2492 

25.0 

4.3274 

50.0 

4. 3941 

75.0 

4.4514 

100.0 

4.5006 

125.0 

4.5430 

150.0 

4.5795 

175.0 

4.6111 

200.0 

4.6384 

225.0 

4.6621 

250.0 

4.6828 

275.0 

4.7008 

300.0 

4.7167 

325.0 

4.7306 

350.0 

4. 7430 

375.0 

4.7539 

400.0 

4.7637 

425.0 

4.7726 

450.0 

4.7805 

475.0 

4.7878 

500.0 

4.7945 

525.0 

4.8006 

550.0 

4.8063 

575.0 

4.8117 

600.0 

4.8167 

625.0 

4.8215 

650.0 

4.8260 

675.0 

4.8303 

700.0 

4.8345 

725.0 

4.8385 

750.0 

4.8425 

775.0 

4.8463 

800.0 

4.8500 

825.0 

4.8537 

850.0 

4.8573 

875.0 

4.8608 

900.  0 

4.8644 

925.0 

4.8678 

950.0 

4.8713 

975.0 

4.8745 

1000.0 

4.8781 

MACH 

RELATIVE 

NUMBER 

HUMIDITY 

0.013322 

0.500 

0.013449 

1.000 

0.013556 

1.000 

0.013649 

1.000 

0.013728 

1.000 

0.013797 

1.000 

0.013857 

1.000 

0.013910 

1.000 

0.013956 

1  .000 

0.013996 

1.000 

0.014032 

1  .000 

0.014063 

1.000 

0.014092 

1.000 

0.014117 

1.000 

0.014140 

1 .000 

0.014161 

1.000 

0.014181 

1.000 

0.014199 

1.000 

0.014215 

1.000 

0.014231 

1.000 

0. 014245 

1.000 

0.014259 

1.000 

0.014272 

1.000 

0.014285 

1.000 

0.014297 

1.000 

0.014309 

1.000 

0.014321 

1.000 

0.014332 

1.000 

0.014343 

1.000 

0.014354 

1.000 

0.014364 

1.000 

0.014375 

1.000 

0.014385 

1 .000 

0.014395 

1.000 

0.014406 

1.000 

0. 014416 

1.000 

0.014426 

1.000 

0.014436 

1.000 

0.014446 

1.000 

0.014455 

1.000 

0.014465 

1  .000 

Case  Jj:  Incompressible  model  -  constrained 


TEMPERATURE  AT  LOWER  ENTRANCE 
PRESSURE  AT  LOWER  ENTRANCE 
PRESSURE  AT  UPPER  ENTRANCE 
RELATIVE  HUMIDITY  AT  LOWER  ENTRANCE 


-20. 000  DEG  C. 
950.000  MB 
925.000  MB 
0.500 
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NUMBER  OF 

DIVISIONS  OF  PASSAGE 

40 

LENGTH  OF 

PASSAGE 

1000. 

000  M 

RADIUS  OF 

PASSAGE 

1. 

500  M 

ANGLE  OF 

PASSAGE 

10. 

000  DEG 

MALL  TEMPERATURE 

10. 

000  DEG  C. 

FRICTION 

FACTOR 

0. 

030 

FRACTION 

OF  HALL  THAT  IS  MOIST 

1. 

000 

DISTANCE 

TEMPERATURE 

PRESSURE 

DENSITY 

(METRES) 

(DEG  CELSIUS) 

(MILLIBARS) 

(KG/CUBIC  M) 

0.0 

-20.00 

949.9 

1.3071 

25.0 

-15.51 

949.2 

1.2833 

50.0 

-11.70 

948.6 

1.2634 

75.0 

-8.46 

947.9 

1.2468 

100.0 

-5.71 

947.2 

1.2328 

125.0 

-3.38 

946.6 

1.2209 

150.0 

-1.40 

945.9 

1.2107 

175.0 

0.28 

945.3 

1.2021 

200.0 

1.71 

944.6 

1.1946 

225.0 

2.92 

944.0 

1.1882 

250.0 

3.95 

943.3 

1.1827 

275.0 

4.82 

942.7 

1.1779 

300.0 

5.57 

942.0 

1.1737 

325.0 

6.20 

941.4 

1.1701 

350.0 

6.73 

940.7 

1 . 1668 

375.0 

7.19 

940.1 

1.1640 

400.  0 

7.57 

939.4 

1.1614 

425.0 

7.90 

938.8 

1.1592 

450.0 

8.18 

938.1 

1.1571 

475.0 

8.41 

937.5 

1.1552 

500.0 

8.61 

936.8 

1 . 1535 

525.0 

8.78 

936.2 

1.1520 

550.0 

8.93 

935.6 

1.1505 

575.0 

9.05 

934.9 

1.1492 

600.0 

9.16 

934.3 

1 .1479 

625.0 

9.24 

933.6 

1.1467 

650.0 

9.32 

933.0 

1.1456 

675.0 

9.38 

932.3 

1.1445 

700.0 

9.44 

931.7 

1.1435 

725.0 

9.48 

931.0 

1.1425 

750.0 

9.52 

930.4 

1.1415 

775.0 

9.55 

929.7 

1.1406 

800. 0 

9.58 

929.1 

1.1397 

825.0 

9.61 

928.5 

1.1388 

85  0.0 

9.63 

927.8 

1.1379 

875.0 

9.64 

927.2 

1.1370 

900.  0 

9.66 

926.5 

1.1362 

925.0 

9.67 

925.9 

1.1353 

950.0 

9.68 

925.2 

1.1345 

975.  0 

9.68 

924.6 

1.1337 

1000.0 

9.70 

925.0 

1. 1329 
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DISTANCE 

VELOCITY 

(METRES) 

(M/SEC) 

0.0 

4.1124 

25.0 

4.1886 

50.0 

4.2544 

75.0 

4.3112 

100.0 

4.3602 

125.0 

4.4027 

150.0 

4.4396 

175.0 

4.4716 

200.0 

4.4995 

225.0 

4.5237 

250.0 

4.5448 

275.0 

4.5633 

300.0 

4.5796 

325.0 

4.5939 

350.0 

4.6066 

375.0 

4.6179 

400.0 

4.6280 

425.0 

4.6371 

450.0 

4.6453 

475.0 

4.6528 

500.0 

4.6597 

525.0 

4.6660 

550.0 

4.6718 

575.0 

4.6773 

600.0 

4.6824 

625.0 

4.6873 

650.0 

4.6919 

675.0 

4.6963 

700.0 

4.7006 

725.0 

4. 7047 

750.0 

4.7087 

775.0 

4.7125 

800.0 

4.7163 

825.0 

4.7200 

850.0 

4.7237 

875.0 

4.  7273 

900.0 

4.7308 

925.0 

4.7344 

950.0 

4.7379 

975.0 

4.7412 

1000.0 

4.7448 

STAGNATION 


PRESSURE 

RELATIVE 

(MILLIBARS) 

HUMIDITY 

950.0 

0.500 

949.3 

1.000 

948.7 

1.000 

948.0 

1  .000 

947.4 

1.000 

946.7 

1.000 

946.0 

1.000 

945.4 

1.000 

944.7 

1.000 

944.1 

1.000 

943.4 

1.000 

942.8 

1.000 

942.1 

1.000 

941.5 

1.000 

940.8 

1.000 

940.2 

1.000 

939.6 

1.000 

938.9 

1.000 

938.3 

1.000 

937.6 

1.000 

937.0 

1.000 

936.3 

1  .000 

935.7 

1.000 

935.0 

1.000 

934.4 

1.000 

933.7 

1.000 

933.1 

1.000 

932.4 

1  .000 

931.8 

1.000 

931  .2 

1  .000 

930.5 

1.000 

929.9 

1.000 

929.2 

1.000 

928.6 

1.000 

927.9 

1.000 

927.3 

1.000 

926.6 

1.000 

926.0 

1.000 

925.4 

1.000 

924.7 

1.000 

925.1 

1.000 

Case  Incompressible  model  -  unconstrained 


TEMPERATURE  AT  LOWEB  ENTEANCE 
PEESSOBE  AT  LOWEB  ENTEANCE 
PRESSURE  AT  UPPER  ENTRANCE 
RELATIVE  HUMIDITY  AT  LCWEE  ENTRANCE 


-20.000  DEG  C. 
950.000  MB 
925.000  MB 
0.500 


.  { 


i 


*> 


V 


' 


100 


NUMBER  OF  DIVISIONS  OF  PASSAGE  40 


LENGTH  OF  PASSAGE 

1000.000 

M 

RADIUS  OF  PASSAGE 

1.500 

M 

ANGLE  OF  PASSAGE 

10.000 

DEG 

WALL  TEMPERATURE 

10.000 

DEG  C. 

FRICTION  FACTOR 

0.030 

FRACTION  OF  WALL  THAT  IS  MOIST 

1.000 

DISTANCE  TEMPERATURE  PRESSURE  DENSITY 

(METRES)  (DEG  CELSIUS)  (MILLIBARS)  (KG/CUBIC  M) 


0.0 

-20.00 

25.0 

-18. 18 

50.0 

-16.23 

75.0 

-14.24 

100.0 

-12.25 

125.0 

-10.31 

150.0 

-8.46 

175.0 

-6.71 

200.0 

-5.07 

225.0 

-3.56 

250.0 

-2.16 

275.0 

-0.88 

300.0 

0.29 

325.0 

0.70 

350.0 

1.71 

375.0 

2.61 

400.0 

3.43 

425.0 

4.  17 

450.0 

4.82 

475.0 

5.41 

500.0 

5.93 

525.0 

6.39 

550.0 

6.80 

575.0 

7.16 

600.0 

7.48 

625.0 

7.76 

650.0 

8.01 

675.0 

8.23 

700.0 

8.42 

725.0 

8.59 

750.0 

8.74 

775.0 

8.87 

800.0 

8.98 

825.0 

9.08 

850.0 

9.17 

875.0 

9.24 

900.0 

9.31 

925.0 

9.37 

950.0 

9.42 

975.  0 

9.46 

1000.0 

9.50 

949.9 

1.3071 

949.2 

1.2961 

948.6 

1.2846 

947.9 

1.2731 

947.2 

1.2620 

946  .6 

1.2514 

945.9 

1.2414 

945.3 

1.2321 

944.6 

1.2234 

944 .0 

1.2155 

943.3 

1.2082 

942.7 

1.2015 

942.0 

1.1954 

941.4 

1.1927 

940.7 

1.1874 

940.1 

1.1826 

939.5 

1.1782 

938.8 

1.1742 

938.2 

1.1706 

937.5 

1.1673 

936.9 

1.1643 

936.2 

1.1615 

935.6 

1.1590 

934.9 

1.1567 

934.3 

1.1546 

933.7 

1.1526 

933.0 

1.1508 

932.4 

1 .1491 

931.7 

1.1475 

931 .1 

1.1460 

930.5 

1. 1446 

929.8 

1.1433 

929.2 

1.1421 

928. S 

1 . 1409 

927.9 

1.1397 

927.3 

1.1386 

926.6 

1.1376 

926 .0 

1.1366 

925.3 

1.1356 

924.7 

1.1346 

925.0 

1.1337 

. 

* 

. -f 
V  .  { 

. 

' 

■ 

.  • 

•  <  < 

I  > 

i 


. 


«  •' 

.  * 


.  • 

► 

. 

■  ♦  ' 

.  * 


r 


'  > 


»  , 


.  I 

« 


.  K 

. 

*• 

t 

, 


•  1 

1  c 

* 


to 


.  . 


, 

.  ) 


. 


DISTANCE 

VELOCITY 

STAGNATION 

PRESSURE 

RELATIVE 

(METRES) 

(M/SEC) 

(MILLIBARS) 

HUMIDITY 

0.0 

4. 0487 

950.0 

0.500 

25.0 

4.0830 

949.3 

2.233 

50.0 

4.  1196 

948.7 

3.1  14 

75.0 

4.1566 

948.0 

3.469 

100.0 

4. 1932 

947.4 

3.499 

125.0 

4.2288 

946.7 

3.355 

150.0 

4.2629 

946.0 

3.131 

175.0 

4.2952 

945.4 

2.880 

200.0 

4.3255 

944.7 

2.634 

225.0 

4.3538 

944.1 

2.406 

250.0 

4.3801 

943.4 

2.202 

275.0 

4.4044 

942.8 

2.023 

300.0 

4.4268 

942.2 

1.868 

325.0 

4.4369 

941.5 

1.739 

350.0 

4.4567 

940.9 

1.719 

375.0 

4.4748 

940.2 

1.623 

400.0 

4.4915 

939.6 

1.540 

425.0 

4.5067 

938.9 

1.468 

450.0 

4.5206 

938.3 

1.406 

475.0 

4.5334 

937.6 

1.352 

500.0 

4.5451 

937.0 

1.305 

525.0 

4.5559 

936.4 

1.265 

550.0 

4.5658 

935.7 

1.230 

575.0 

4.5750 

935.1 

1.199 

600.0 

4.5834 

934.4 

1.  173 

625.0 

4.5912 

933.8 

1.150 

650.0 

4.5985 

933.1 

1.129 

675.0 

4.6053 

932.5 

1.112 

700.0 

4.6116 

931  .9 

1.096 

725.0 

4.6176 

931 .2 

1.083 

750.0 

4.6232 

930.6 

1.071 

775.0 

4.6286 

929.9 

1.061 

800.0 

4.6336 

929.3 

1.052 

825.0 

4.6384 

928.7 

1.043 

850.0 

4.6431 

928.0 

1.036 

875.0 

4.6475 

927.4 

1.030 

900.0 

4.6518 

926.7 

1.025 

925.0 

4.6560 

926.1 

1.020 

950.0 

4.6600 

925.5 

1.015 

975.0 

4.6640 

924.8 

1.012 

1000.0 

4.6678 

925.1 

1.008 

Case  2i  Compressible  model 

TEMPERATURE  AT  LOWER  ENTRANCE 

-10 

.000  BEG  C 

PRESSURE  AT 

LOWER  ENTRANCE 

950 

.000  MB 

PRESSURE  AT 

UPPER  ENTRANCE 

925 

.000  MB 

RELATIVE  HUMIDITY  AT  LOWER 

ENTRANCE 

0 

.500 
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NUMBER  OF  DIVISIONS  OF  PASSAGE  40 


LENGTH  OF  PASSAGE 

1000.000 

M 

RADIUS  OF  PASSAGE 

1.500 

M 

ANGLE  OF  PASSAGE 

10.000 

DEG 

WALL  TEMPERATURE 

10.000 

DEG  C. 

FRICTION  FACTOR 

0.030 

FRACTION  OF  WALL  THAT  IS  MOIST 

1.000 

DISTANCE 

TEMPERATURE 

(METRES) 

(DEG  CELSIUS) 

0.0 

-9.99 

25.0 

-7.02 

50.0 

-4.50 

75.0 

-2.36 

100.0 

-0.54 

125.0 

1.01 

150.0 

2.32 

175.0 

3.44 

200.0 

4.39 

225.0 

5.19 

250.0 

5.87 

275.0 

6.46 

300.0 

6.95 

325.0 

7.37 

350.0 

7.73 

375.0 

8.03 

400.0 

8.29 

425.0 

8.51 

450.0 

8.6  9 

475.0 

8.85 

500.0 

8.98 

525.0 

9.10 

550.0 

9.19 

575.0 

9.28 

600.0 

9.35 

625.0 

9.40 

650.0 

9.46 

675.0 

9.50 

700.0 

9.53 

725.0 

9.57 

750.0 

9.59 

775.0 

9.61 

800.0 

9.63 

825.0 

9.65 

850.  0 

9.66 

875.0 

9.67 

900.0 

9.67 

925.0 

9.69 

950.0 

9.69 

975.0 

9.71 

1000.0 

9.71 

PRESSURE 

DENSITY 

(MILLIBARS) 

(KG/CUBIC  M) 

950.0 

1.2576 

949.4 

1.2427 

948.8 

1.2301 

948.1 

1.2196 

947.5 

1.2106 

946.9 

1.2029 

946.2 

1.1963 

945.6 

1.1907 

945.0 

1.1858 

944.4 

1.1816 

943.7 

1.1779 

943.1 

1.1746 

942.5 

1.1717 

941 .8 

1.1692 

941.2 

1 .1669 

940.6 

1.1649 

940.0 

1.1630 

939.3 

1.1613 

938.7 

1.1598 

938.1 

1.1583 

937.4 

1.1570 

936.8 

1.1558 

936.2 

1.1546 

935.6 

1.1535 

935.0 

1.1524 

934.3 

1.1514 

933.7 

1.1504 

933.1 

1.1495 

932.5 

1.1486 

931.8 

1.1477 

931.2 

1.1468 

930.6 

1.1460 

930.0 

1.1451 

929.3 

1.1443 

928.7 

1 .1435 

928.1 

1.1427 

927  .5 

1.1419 

926.9 

1.1411 

926.2 

1 .1403 

925.6 

1.1395 

925  .0 

1.1387 
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STAGNATION 


DISTANCE 

TEMPERATURE 

(METRES) 

(DEG  CELSIUS) 

0.0 

-9.98 

25.0 

-7.01 

50.0 

-4.49 

75.0 

-2.35 

100.0 

-0.53 

125.  0 

1.02 

150.0 

2.33 

175.0 

3.45 

200.0 

4.40 

225.0 

5.20 

250.0 

5.89 

275.0 

6.47 

300.0 

6.96 

325.0 

7.38 

350.0 

7.74 

375.0 

8.04 

400.0 

8.30 

425.0 

8.52 

450.0 

8.70 

475.0 

8.86 

500.0 

8.99 

525.  0 

9.11 

550.0 

9.20 

575.0 

9.29 

600.0 

9.36 

625.  0 

9.42 

650.0 

9.47 

675.0 

9.51 

700.0 

9.55 

725.0 

9.58 

750.0 

9.60 

775.0 

9.63 

800.0 

9.64 

825.0 

9.66 

850.0 

9.67 

875.0 

9.69 

900.0 

9.69 

925.0 

9.70 

950.0 

9.70 

975.  0 

9.72 

1000. 0 

9.72 

DISTANCE 

VELOCITY 

(METRES) 

(M/SEC) 

0.0 

4.4200 

25.0 

4.4730 

50.0 

4.5185 

75.0 

4.5577 

100.0 

4.5915 

STAGNATION 

STAGNATION 

PRESSURE 

DENSITY 

(MILLIBARS) 

(KG/CUBIC  M) 

950.1 

1.2577 

949.5 

1.2428 

948.9 

1.2303 

948.3 

1.2197 

947.6 

1.2107 

947.0 

1.2030 

946.4 

1.1965 

945.7 

1.1908 

945.1 

1.1859 

944.5 

1.1817 

943.9 

1.1780 

943.2 

1.1747 

942.6 

1.1718 

942.0 

1.1693 

941 .3 

1.1670 

940.7 

1.1650 

940.1 

1.1631 

939.5 

1.1614 

938.8 

1.1599 

938.2 

1.1585 

937.6 

1.1571 

937.0 

1.1559 

936.3 

1.1547 

935.7 

1. 1536 

935.1 

1.1526 

934.5 

1.1515 

933.8 

1.1506 

933.2 

1 . 1496 

932.6 

1.1487 

932.0 

1.1478 

931.3 

1.1469 

930.7 

1. 1461 

930.1 

1.1452 

929.5 

1. 1444 

928.9 

1.1436 

928.2 

1. 1428 

927.6 

1.1420 

927.0 

1.1412 

926.4 

1.1404 

925.8 

1.1396 

925.1 

1.1388 

MACH 

RELATIVE 

NUMBER 

HUMIDITY 

0.013591 

0.499 

0.013677 

1.000 

0.013751 

1  .000 

0.013815 

1.000 

0.013871 

1.000 
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125.0 

4.6208 

0.013919 

1.000 

150.0 

4.6462 

0.013962 

1.000 

175.0 

4.6682 

0.014000 

1 .000 

200.0 

4.6875 

0.014034 

1.000 

225.0 

4. 7043 

0.014063 

1.000 

250.0 

4.7191 

0.014090 

1.000 

275.0 

4.7322 

C. 0141 14 

1.000 

300.0 

4.7438 

0.014136 

1.000 

325.0 

4.7541 

0.014157 

1.000 

350.0 

4.7634 

0.014175 

1.000 

375.0 

4.7718 

0.014192 

1.000 

400.0 

4.7794 

0.014208 

1.000 

425.0 

4.7863 

0.014223 

1.000 

450.0 

4.7927 

0.014238 

1.000 

475.0 

4.  7986 

0.014251 

1.000 

500.0 

4.8041 

0.01 4264 

1.000 

525.0 

4.8093 

0.014276 

1.000 

550.0 

4.8141 

0.014288 

1.000 

575.0 

4.8188 

0. 014300 

1  .000 

600.0 

4.8232 

0.014311 

1.000 

625.0 

4.8274 

0.014323 

1.000 

650.0 

4.8315 

0.014333 

1.000 

675.0 

4.8355 

0.014344 

1 .000 

7C0.0 

4.8393 

0.014355 

1.000 

725.0 

4.8431 

0.014365 

1.000 

750.0 

4.8468 

0.014375 

1.000 

775.0 

*  4.8504 

0.014385 

1 .000 

800.0 

4.8540 

0.014396 

1.000 

825.0 

4.8575 

0.014406 

1.000 

850.0 

4.8610 

0.014416 

1.000 

875.0 

4.8645 

0.014426 

1.000 

900.0 

4.8677 

0.014435 

1.000 

925.0 

4.8713 

0.014445 

1.000 

950-0 

4.8746 

0.014455 

1.000 

975.0 

4.8781 

0.014465 

1.000 

1C00.0 

4.8813 

0.014475 

1.000 

Case  2: 

Inccmpressible 

model  -  constrained 

TEMPERATURE  AT  LOWER  ENTRANCE 
PRESSURE  AT  LOWER  ENTRANCE 
PRESSURE  AT  UPPER  ENTRANCE 
RELATIVE  HUMIDITY  AT  LOWER  ENTRANCE 
NUMBER  OE  DIVISIONS  OF  PASSAGE 


-10.000  DEG  C 
950.000  MB 
925.000  MB 
0.500 
40 


LENGTH  OF  PASSAGE 
RADIUS  OF  PASSAGE 
ANGLE  OF  PASSAGE 
WALL  TEMPERATURE 
FRICTION  FACTOR 

FRACTION  OF  WALL  THAT  IS  MOIST 


1000.000  M 
1.500  M 
10.000  DEG 
10.000  DEG  C 
0.030 
1.000 


r 


, 


*  ' 
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DISTANCE 

TEMPERATURE 

(METRES) 

(DEG  CELSIUS) 

0.0 

-10.00 

25.0 

-7.02 

50.0 

-4.50 

75.0 

-2.35 

100.0 

-0.53 

125.0 

1.02 

150.0 

2.34 

175.0 

3.45 

200.0 

4.40 

225.0 

5.21 

250.0 

5.89 

275.0 

6.47 

300.0 

6.96 

325.0 

7.38 

350.0 

7.74 

375.0 

8.04 

400.0 

8.30 

425.0 

8.51 

450.0 

8.70 

475.0 

8.86 

500.0 

8.99 

525.0 

9.10 

550.0 

9.20 

575.0 

9.28 

600.0 

9.35 

625.0 

9.41 

650.0 

9.46 

675.0 

9.50 

7C0.0 

9.54 

725.0 

9.57 

750.0 

9.59 

775.0 

9.62 

800.0 

9.64 

825.0 

9.65 

850.0 

9.66 

875.0 

9.68 

900.0 

9.68 

925.0 

9.69 

950.0 

9.69 

975.0 

9.71 

1C00.0 

9.71 

DISTANCE 

VELOCITY 

(METRES) 

(M/SEC) 

0.0 

4.3200 

25.0 

4.3736 

50.0 

4.4195 

75.0 

4.4593 

100.0 

4.4939 

PRESSURE 

DENSITY 

(MILLIBARS) 

(KG/CUBIC  M) 

949.9 

1.2570 

949.2 

1.2416 

948.6 

1.2287 

947.9 

1.2177 

947.3 

1.2084 

946.6 

1.2003 

946.0 

1.1934 

945.3 

1.1875 

944.7 

1. 1824 

944.1 

1.1779 

943.4 

1.1740 

942.8 

1.1705 

942.1 

1.1675 

941.5 

1.1648 

940.8 

1.1624 

940.2 

1.1603 

939.6 

1.1583 

938.9 

1.1565 

938.3 

1.1549 

937.6 

1.1534 

937.0 

1.1520 

936.4 

1.1507 

935.7 

1.1495 

935.1 

1.1484 

934.4 

1.1473 

933.8 

1.1462 

933.2 

1.1452 

932.5 

1.1442 

931.9 

1.1433 

931.2 

1.1424 

930.6 

1.1415 

930.0 

1.1406 

929.3 

1.1397 

928.7 

1.1389 

928.0 

1.1380 

927.4 

1.1372 

926.8 

1.1364 

926.1 

1.1356 

925.5 

1.1348 

924.8 

1.1339 

925.0 

1.1331 

STAGNATION 

PRESSURE 

RELATIVE 

(MILLIBARS) 

HUMIDITY 

950.0 

0.500 

949.3 

1.000 

948.7 

1.000 

948.0 

1.000 

947.4 

1.000 

** 
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125.0 

4.5240 

150.0 

4.5501 

175.0 

4.5728 

200.0 

4.5927 

225.0 

4.6101 

250.0 

4.6255 

275.0 

4.6391 

300.0 

4.6511 

325.0 

4.6619 

350.0 

4.6715 

375.0 

4.6802 

400.0 

4.6880 

425.0 

4.6952 

450.0 

4.7018 

475.0 

4.7079 

500.0 

4. 7136 

525.0 

4.7189 

550.0 

4.7239 

575.0 

4.7286 

600.0 

4.7332 

625.0 

4.7375 

650.0 

4.7417 

675.0 

4.7457 

700.0 

4.7496 

725.0 

4.7535 

750.0 

4.7572 

775.0 

4.7609 

800.0 

4.7645 

825.0 

4.7681 

850.0 

4.7716 

875.0 

4.7751 

900.0 

4.7784 

925.0 

4.7820 

950.0 

4.7854 

975.0 

4.7889 

1000.0 

4.7922 

946.8 

1.000 

946.1 

1.000 

945.5 

1.000 

944.8 

1.000 

944.2 

1.000 

943.5 

1.000 

942.9 

1.000 

942.3 

1.000 

941.6 

1.000 

941.0 

1 .000 

940.3 

1.000 

939.7 

1.000 

939.0 

1.000 

938.4 

1.000 

937.8 

1.000 

937.1 

1.000 

936.5 

1.000 

935.8 

1.000 

935.2 

1.000 

934.6 

1.000 

933.9 

1.000 

933.3 

1.000 

932.6 

1.000 

932.0 

1  .000 

931.4 

1.000 

930.7 

1.000 

930.1 

1.000 

929.4 

1  .000 

928.8 

1.000 

928.2 

1.000 

927.5 

1.000 

926.9 

1.000 

926.3 

1.000 

925.6 

1.000 

925.0 

1.000 

925.1 

1.000 

Case  2:  Incompressible  model  -  unconstrained 


TEMPERATURE  AT  LOWER  ENTRANCE 
PRESSURE  AT  LOWER  ENTRANCE 
PRESSURE  AT  UPPER  ENTRANCE 
RELATIVE  HUMIDITY  AT  LOWER  ENTRANCE 
NUMBER  OF  DIVISIONS  OF  PASSAGE 


-10.000  DEG  C 
950.000  MB 
925.000  MB 
0.500 
40 


LENGTH  OF  PASSAGE 
RADIUS  OF  PASSAGE 
ANGLE  OF  PASSAGE 
WALL  TEMPERATURE 
FRICTION  FACTOR 

FRACTION  OF  WALL  THAT  IS  MOIST 


1000.000  M 
1.500  M 
10.000  DEG 
10.000  DEG  C 
0.030 
1.000 
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DISTANCE 

TEMPERATURE 

(METRES) 

(DEG  CELSIUS) 

0.0 

-10.00 

25.0 

-9.52 

50.0 

-8.73 

75.0 

-7.74 

100.0 

-6.63 

125.0 

-5.45 

150.0 

-4.26 

175.0 

-3.08 

200.0 

-1.94 

225.0 

-0.85 

250.0 

0.18 

275.0 

0.43 

300.0 

1.37 

325.0 

2.24 

350.0 

3.04 

375.0 

3.76 

400.  0 

4.42 

425.0 

5.02 

450.0 

5.56 

475.0 

6.04 

500.0 

6.47 

525.0 

6.86 

550.0 

7.21 

575.0 

7.51 

600.0 

7.78 

625.0 

8.02 

650.0 

8.24 

675.0 

8.42 

700.0 

8.59 

725.0 

8.73 

750.0 

8.86 

775.0 

8.97 

800.  0 

9.07 

825.0 

9.16 

850.0 

9.23 

875.0 

9.30 

900.0 

9.36 

925.0 

9.41 

950.0 

9.45 

975.0 

9.49 

1000.0 

9.52 

DISTANCE 

VELOCITY 

(METRES) 

(M/SEC) 

0.0 

4.2506 

25.0 

4.2636 

50.0 

4.2817 

75.0 

4.3027 

100.0 

4.3254 

PRESSURE 

DENSJTY 

(MILLIBARS) 

(KG/CUBIC  M) 

949.9 

1.2570 

949.2 

1.2532 

948.6 

1.2479 

947.9 

1.2418 

947.3 

1.2353 

946.6 

1.2286 

946.0 

1.2219 

945.4 

1.2155 

944.7 

1.2093 

944.1 

1.2034 

943.4 

1.1979 

942.8 

1.1958 

942.1 

1.1908 

941.5 

1.1861 

940.9 

1.1818 

940.2 

1.1778 

939.6 

1.1742 

939.0 

1.1708 

938.3 

1.1677 

937.7 

1.1648 

937.0 

1.1622 

936.4 

1.1598 

935.8 

1. 1576 

935.1 

1.1555 

934.5 

1.1536 

933.9 

1.1518 

933.2 

1 . 1501 

932.6 

1.1486 

932.0 

1.1471 

931.3 

1.1457 

930.7 

1.  1444 

930.0 

1.1432 

929.4 

1.1420 

928.8 

1.1409 

928.1 

1.1398 

927.5 

1.1387 

926.9 

1.1377 

926.2 

1.1367 

925.6 

1. 1358 

925.0 

1. 1348 

925.0 

1.1339 

STAGNATION 

PRESSURE 

RELATIVE 

(MILLIBARS) 

HUMIDITY 

950.0 

0.500 

949.4 

1.141 

948.7 

1.616 

948.1 

1.920 

947.4 

2.081 

' 
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125.0 

4.3489 

150.0 

4.3726 

175.0 

4.3958 

200.0 

4.4183 

225.0 

4.4399 

250.0 

4.4603 

275.0 

4.4681 

300.0 

4.4869 

325.  0 

4.5046 

350.0 

4.5210 

375.0 

4.5363 

400.0 

4.5505 

425.0 

4.5636 

450.0 

4.5757 

475.0 

4.5869 

500.0 

4.5973 

525.0 

4.6069 

550.0 

4.6158 

575.0 

4.6240 

600.0 

4.6317  - 

625.0 

4.6389 

650.0 

4.6456 

675.0 

4.6519 

700.0 

4.6578 

725.0 

4.6635 

750.0 

4.6688 

775.0 

4.6738 

800.0 

4.6787 

825.0 

4.6833 

850.0 

4.6878 

875.0 

4.6921 

900.0 

4.6963 

925.0 

4.7003 

950.0 

4.7043 

975.0 

4.7081 

1000.0 

4.7119 

946.8 

2.137 

946.1 

2.120 

945.5 

2.058 

944.8 

1.973 

944  .2 

1.877 

943.5 

1.779 

942.9 

1.688 

942.3 

1.699 

941  .6 

1.621 

941.0 

1.549 

940.3 

1.484 

939.7 

1.426 

939.1 

1.374 

938.4 

1.328 

937.8 

1.288 

937.2 

1.252 

936.5 

1.220 

935.9 

1.193 

935.3 

1.168 

934.6 

1.146 

934.0 

1.128 

933.3 

1.111 

932.7 

1.096 

932.1 

1.083 

931  .4 

1.072 

930.8 

1.062 

930.2 

1.053 

929.5 

1.045 

928.9 

1.038 

928.3 

1.031 

927.6 

1.026 

927.0 

1.021 

926.4 

1.017 

925.7 

1.013 

925.1 

1.010 

925.1 

1.006 

Case  3 :  Compressible  model 


TEMPERATURE  AT  LOWER  ENTRANCE 
PRESSURE  AT  LOWER  ENTRANCE 
PRESSURE  AT  UPPER  ENTRANCE 
RELATIVE  HUMIDITY  AT  LCWER  ENTRANCE 
NUMBER  OE  DIVISIONS  OF  PASSAGE 


-.0  DEG  C 
950.000  MB 
925.000  MB 
0.500 
40 


LENGTH  OF  PASSAGE 
RADIUS  OF  PASSAGE 
ANGLE  OF  PASSAGE 
WALL  TEMPERATURE 
FRICTION  FACTOR 

FRACTION  OF  WALL  THAT  IS  MOIST 


1000.000  M 
1.500  M 
10.000  DEG 
10.000  DEG  C 
0.030 
1.000 
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DISTANCE 

TEMPEEATUBE 

(METEES) 

(DEG  CELSIUS) 

0.0 

0.01 

25.0 

-0.89 

50.0 

-0.84 

75.0 

0.75 

100.0 

2.11 

125.0 

3.26 

150.0 

4.23 

175.0 

5.06 

200.0 

5.76 

225.0 

6.36 

250.0 

6.87 

275.0 

7.30 

300.0 

7.67 

325.0 

7.98 

350.0 

8.25 

375.0 

8.47 

400.0 

8.66 

425.0 

8.82 

450.0 

8.96 

475.  0 

9.08 

500.0 

9.18 

525.0 

9.26 

550.0 

9.33 

575.0 

9.40 

600.0 

9.45 

625.0 

9.49 

650.0 

9.53 

675.  0 

9.56 

700.0 

9.59 

725.0 

9.61 

750.0 

9.63 

775.0 

9.65 

800.0 

9.66 

825.0 

9.67 

850.0 

9.68 

875.  0 

9.68 

900.0 

9.70 

925.0 

9.70 

950.0 

9.71 

975.0 

9.71 

1 000.0 

9.71 

STAGNATION 

DISTANCE 

TEMPEEATUEE 

(METEES) 

(DEG  CELSIUS) 

0.0 

0.02 

25.0 

-0.88 

50.0 

-0.83 

75.0 

0.76 

100.0 

2.12 

PEESSUEE 

DENSITY 

(MILLIB AES) 

(KG/CUBIC  M) 

950.0 

1.2114 

949.4 

1.2146 

948.7 

1.2134 

948.1 

1.2056 

947.5 

1.1988 

946.8 

1.1930 

946.2 

1.1880 

945.6 

1.1836 

944.9 

1.1798 

944  .3 

1.1765 

943.7 

1. 1736 

943.1 

1.1710 

942.4 

1.1687 

941  .8 

1.1666 

941.2 

1.1647 

940.6 

1.1630 

939.9 

1.1614 

939.3 

1.1600 

938.7 

1.1586 

938.1 

1.1574 

937.4 

1.1562 

936.8 

1.1551 

936.2 

1.1540 

935.6 

1.1530 

934.9 

1.1520 

934.3 

1.1510 

933.7 

1.1501 

933.1 

1.1492 

932.4 

1.1484 

931  .8 

1.1475 

931.2 

1.1467 

930.6 

1.1458 

930.0 

1.1450 

929.3 

1.1442 

928.7 

1. 1434 

928.1 

1.1426 

927.5 

1.1418 

926  .9 

1.1410 

926.2 

1.1402 

925.6 

1.1395 

925.0 

1.1387 

STAGNATION 

STAGNATION 

PEESSUEE 

DENSITY 

(MIL LIE AES) 

(KG/CUBIC  M) 

950.1 

1.2115 

949.5 

1.2147 

948.9 

1.2136 

948.2 

1.2057 

947.6 

1.1989 
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125.0 

3.27 

150.0 

4.24 

175.0 

5.07 

200.0 

5.78 

225.0 

6.37 

250.0 

6.88 

275.0 

7.31 

300.0 

7.68 

325.0 

7.99 

350.0 

8.26 

375.0 

8.48 

400.0 

8.67 

425.0 

8.84 

450.0 

8.97 

475.0 

9.09 

500.0 

9.19 

525.0 

9.27 

550.0 

9.35 

575.0 

9.41 

600.0 

9.46 

625.0 

9.50 

650.0 

9.54 

675.0 

9.57 

700.  0 

9.60 

725.0 

9.62 

750.0 

9.64 

775.0 

9.66 

800.0 

9.67 

825.0 

9.68 

850.  0 

9.69 

875.0 

9.69 

900.0 

9.71 

925.0 

9.71 

950.  0 

9.72 

975.0 

9.72 

1000.0 

9.72 

DISTANCE 

VELOCITY 

(METRES) 

( M/SEC) 

0.0 

4.5701 

25.0 

4.5583 

50.0 

4.5626 

75.0 

4.5923 

100.0 

4.6182 

125.0 

4.6407 

150.0 

4.6603 

175.0 

4.6774 

200.0 

4.6925 

225.0 

4.7057 

250.0 

4.7175 

275.0 

4.7280 

300.0 

4.7374 

325.0 

4.7458 

947.0 

1.1931 

946.3 

1.1881 

945.7 

1.1838 

945.1 

1.1800 

944.4 

1.1766 

943.8 

1.1737 

943.2 

1.1711 

942.6 

1.1688 

941  .9 

1.1667 

941.3 

1.1648 

940.7 

1.1631 

940.1 

1.1615 

939.4 

1.1601 

938.8 

1.1587 

938.2 

1.1575 

937.6 

1.1563 

936.9 

1.1552 

936.3 

1.1541 

935.7 

1.1531 

935.1 

1.1521 

934.4 

1.1512 

933.8 

1.1502 

933.2 

1.1493 

932.6 

1.1485 

932.0 

1.1476 

931.3 

1.1468 

930.7 

1.1459 

930.1 

1.1451 

929.5 

1. 1443 

928.9 

1.1435 

928.2 

1.1427 

927.6 

1.1419 

927.0 

1.1411 

926.4 

1.1403 

925.8 

1.1396 

925.1 

1.1388 

MACH 

RELATIVE 

NUMBER 

HUMIDITY 

0.013793 

0.500 

0.013779 

0.727 

0.013791 

0.993 

0.013840 

1.000 

0.013884 

1.000 

0.013922 

1.000 

0.013956 

1.000 

0.013986 

1.000 

0.014013 

1.000 

0.014038 

1.000 

0.014060 

1.000 

0.014080 

1.000 

0.014099 

1.000 

0.014116 

1.000 
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350.0 

4.7535 

0.014133 

1  .000 

375.0 

4.7605 

0.014148 

1.000 

400.0 

4. 7669 

0.014162 

1.000 

425.0 

4.7729 

0.014175 

1.000 

450.0 

4.7784 

0.014188 

1.000 

475.0 

4.7836 

0.014201 

1.000 

500.0 

4.7885 

0.014213 

1 .000 

525.0 

4.7931 

0.014224 

1.000 

550.0 

4.7976 

0.014236 

1.000 

575.0 

4.8018 

0.014247 

1.000 

600.0 

4.8059 

0.014258 

1.000 

625.0 

4.8099 

0.014268 

1.000 

650.0 

4.8137 

0.014279 

1.000 

675.0 

4.8175 

0.014289 

1.000 

700.0 

4.8212 

0.014299 

1.000 

725.0 

4.8248 

0.014309 

1.000 

750.0 

4.8283 

0.014319 

1  .000 

775.0 

4.8318 

0.014329 

1.000 

800.0 

4.8353 

0.014339 

1.000 

825.0 

4.8387 

0.014349 

1.000 

850.0 

4.8421 

0.014359 

1.000 

875.0 

4.8454 

0.014369 

1.000 

900.0 

4.8489 

0.014379 

1.000 

925.0 

4.8521 

0.014388 

1.000 

950.0 

4.8556 

0.014398 

1.000 

975.0 

4.8588 

0.014408 

1.000 

1000.0 

4.8621 

0.014417 

1.000 

Case  3:  Incompressible 

model  -  constrained 

TEMPERATURE 

AT  LOWER  ENTRANCE 

-.0  DEG  C. 

PRESSURE 

AT 

LOWER  ENTRANCE 

950.000  MB 

PRESSURE 

AT 

UPPER  ENTRANCE 

925.000  MB 

RELATIVE 

HUMIDITY  AT  LOWER 

ENTRANCE 

0.500 

NUMBER  OP 

DIVISIONS  OF 

PASSAGE 

40 

LENGTH  OF 

PASSAGE 

1000.000  M 

RADIUS  OF 

PASSAGE 

1.500  M 

ANGLE  OF 

PASSAGE 

10.000  DEG 

WALL  TEMPERATURE 

10.000  DEG  C. 

FRICTION 

FACTOR 

0.030 

FRACTION 

OF 

WALL  THAT 

IS  MOIST 

1.000 

DISTANCE 

TEMPERATURE 

PRESSURE 

DENSITY 

(METRES) 

(DEG  CELSIUS) 

(MILLIBARS) 

(KG/CUBIC  M) 

0.0 

0.0 

949.9 

1.2101 

25.0 

-0.90 

949.2 

1.2128 

50.0 

-0.84 

948.6 

1.2112 

75.0 

0.75 

948.0 

1.2033 

100.0 

2.11 

947.3 

1.1962 

125.  0 

3.26 

946.7 

1.1901 
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150.0 

4.23 

946.0 

1.1848 

175.0 

5.06 

945.4 

1.1802 

200.  0 

5.77 

944.8 

1.  1762 

225.0 

6.37 

944.1 

1.1727 

250.0 

6.88 

943.5 

1.1696 

275.0 

7.31 

942.9 

1.1669 

300.0 

7.67 

942.2 

1.  1644 

325.0 

7.99 

941  .6 

1.1622 

350.0 

8.25 

941.0 

1.1603 

375.0 

8.47 

940.3 

1.1585 

400.0 

8.67 

939.7 

1. 1568 

425.0 

8.83 

939.1 

1.1553 

450.  0 

8.97 

938.4 

1.1539 

475.0 

9.08 

937.8 

1.1526 

500.0 

9.18 

937.2 

1.1514 

525.0 

9.27 

936.5 

1.1502 

550.0 

9.34 

935.9 

1.1491 

575.0 

9.40 

935.3 

1.1481 

600.0 

9.45 

934  .6 

1.1471 

625.0 

9.49 

934.0 

1.1461 

650.  0 

9.53 

933.4 

1.1452 

675.0 

9.56 

932.7 

1.1442 

700.0 

9.59 

932.1 

1.  1434 

725.0 

9.61 

931.5 

1.1425 

750.0 

9.63 

930.8 

1.1416 

775.0 

9.65 

930.2 

1.1408 

800.0 

9.66 

929.6 

1.1399 

825.0 

9.67 

929.0 

1.1391 

850.0 

9.68 

928.3 

1.1383 

875.0 

9.68 

927.7 

1.1375 

900.0 

9.70 

927.1 

1.1367 

925.0 

9.70 

926.4 

1.1359 

950.0 

9.71 

925.8 

1.1351 

975.0 

9.71 

925.2 

1.1343 

1000.0 

9.71 

925.0 

1. 1335 

STAGNATION 

DISTANCE 

VELOCITY 

PRESSURE 

RELATIVE 

(METRES) 

(H/SEC) 

(MILLIBARS) 

HUMIDITY 

0.0 

4.5100 

950.0 

0.500 

25.0 

4.5000 

949.4 

0.728 

50.0 

4.5062 

948.7 

0.995 

75.0 

4.5357 

948.1 

1.000 

100.0 

4.5626 

947.4 

1.000 

125.0 

4.5860 

946.8 

1.000 

150.0 

4.6064 

946.2 

1.000 

175.0 

4.6243 

945.5 

1.000 

200.0 

4.6400 

944.9 

1  .000 

225.0 

4.6539 

944.3 

1.000 

250.0 

4.6662 

943.6 

1.000 

275.0 

4.6772 

943.0 

1.000 

300.0 

4.6870 

942.4 

1.000 

325.0 

4.6959 

941.7 

1.000 
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350.0 

4.7039 

941.1 

1.000 

375.0 

4.7111 

940.5 

1.000 

400.0 

4.7178 

939.8 

1.000 

425.0 

4.7240 

939.2 

1 .000 

450.0 

4.7297 

938.6 

1.000 

475.0 

4.7351 

937.9 

1.000 

500.0 

4.7401 

937.3 

1.000 

525.0 

4.7449 

936.7 

1.000 

550.0 

4.7494 

936.0 

1.000 

575.0 

4.7537 

935.4 

1.000 

600.0 

4.7579 

934.8 

1.000 

625.0 

4.7620 

934.1 

1.000 

650.0 

4.7659 

933.5 

1.000 

675.0 

4.7697 

932.9 

1 .000 

700.0 

4.7734 

932.2 

1.000 

725.0 

4.7771 

931.6 

1.000 

750.0 

4.7807 

931.0 

1.000 

775.0 

4. 7842 

930.3 

1.000 

800.0 

4.7877 

929.7 

1.000 

825.0 

4.7912 

929.1 

1.000 

850.0 

4.7947 

928.5 

1.000 

875.0 

4.7980 

927.8 

1.000 

900.0 

4.8015 

927.2 

1.000 

925.0 

4.8048 

926.6 

1.000 

950.0 

4.8083 

925.9 

1.000 

975.0 

4.8115 

925.3 

1.000 

1000.0 

4.8148 

925.1 

1.000 

Case  3j_  Incompressible  model  -  unconstrained 


TEMPERATURE  AT  LOWER  ENTRANCE 
PRESSURE  AT  LOWER  ENTRANCE 
PRESSURE  AT  UPPER  ENTRANCE 
RELATIVE  HUMIDITY  AT  LOWER  ENTRANCE 
NUMBER  OF  DIVISIONS  OF  PASSAGE 


-.0  DEG  C. 
950.000  MB 
925.000  MB 
0.500 
40 


LENGTH  OF 

PASSAGE 

1COO.OOO  M 

RADIUS  OF  PASSAGE 

ANGLE  OF  PASSAGE 

WALL  TEMPERATURE 

FRICTION  FACTOR 

FRACTION  OF  WALL  THAT  IS  MOIST 

1.500  M 

10.000  DEG 
10.000  DEG  C. 
0.030 

1.000 

DISTANCE 

TEMPERATURE 

PRESSURE 

DENSITY 

(METRES) 

(DEG  CELSIUS) 

(MILLIBARS) 

(KG/CUBIC  M) 

0.0 

0.0 

949.9 

1.2101 

25.0 

-0.90 

949.2 

1.2128 

50.0 

-0.84 

948.6 

1.2112 

75.0 

-0.76 

948.0 

1.2095 

100.0 

-0.46 

947.3 

1.2069 

125.0 

-0.01 

946.7 

1.2038 
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150.0 

0.54 

175.0 

0.41 

200.0 

1.07 

225.0 

1.73 

250.0 

2.40 

275.0 

3.04 

300.0 

3.66 

325.0 

4.24 

350.0 

4.79 

375.0 

5.30 

400.0 

5.76 

425.0 

6.19 

450.0 

6.58 

475.0 

6.93 

500.0 

7.25 

525.0 

7.53 

550.0 

7.79 

575.0 

8.02 

600.0 

8.22 

625.0 

8.41 

650.0 

8.57 

675.0 

8.71 

700.0 

8.84 

725.0 

8.95 

750.0 

9.05 

775.0 

9.13 

800.0 

9.21 

825.0 

9.28 

850.0 

9.34 

875.0 

9.39 

900.0 

9.43 

925.0 

9.47 

950.0 

9.51 

975.0 

9.54 

1000.0 

9.57 

DISTANCE 

VELOCITY 

(METRES) 

(M/SEC) 

0.0 

4.4700 

25.0 

4.4601 

50.0 

4.4662 

75.0 

4.4724 

100.0 

4.4819 

125.0 

4.4936 

150.0 

4.5066 

175.  0 

4.5086 

200.0 

4.5232 

225.0 

4.5379 

250.0 

4.5525 

275.0 

4.5667 

300.0 

4.5804 

325.0 

4.5935 

946.1 

1  .2003 

945.4 

1.1998 

944.8 

1.1959 

944.2 

1.1920 

943.5 

1.1882 

942.9 

1.1845 

942.3 

1.1810 

941.6 

1.1776 

941.0 

1.1744 

940.4 

1.1714 

939.7 

1. 1686 

939.1 

1.1660 

938.5 

1.1636 

937.8 

1.1613 

937.2 

1.1592 

936.6 

1.1572 

935.9 

1.1553 

935.3 

1.1536 

934.7 

1.1520 

934  .0 

1.1504 

933.4 

1.1490 

932.8 

1.1476 

932.2 

1.1463 

931 .5 

1.1451 

930.9 

1.1439 

930.3 

1.1428 

929.6 

1.1417 

929.0 

1.1406 

928.4 

1.1396 

927.7 

1.1387 

927.1 

1.1377 

926.5 

1.1368 

925.9 

1.1359 

925.2 

1.1350 

925.0 

1.1341 

STAGNATION 

PRESSURE 

RELATIVE 

(MILLIBARS) 

HUMIDITY 

950.0 

0.500 

949.4 

0.728 

948.7 

0.994 

948.1 

1.165 

947.5 

1.306 

946.8 

1.396 

946.2 

1.446 

945.5 

1.471 

944.9 

1.554 

944  .3 

1.538 

943.6 

1.511 

943.0 

1.477 

942.4 

1.440 

941.7 

1.401 
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350.0 

4.6059 

941 .1 

1.364 

375.0 

4.6177 

940.5 

1.328 

400.0 

4.6288 

939.8 

1.294 

425.0 

4.6391 

939.2 

1.262 

450.0 

4.6489 

938.6 

1.234 

475.0 

4.6580 

938.0 

1.207 

500.0 

4.6665 

937.3 

1.184 

525.0 

4.6745 

936.7 

1.162 

550.0 

4.6820 

936 .1 

1.143 

575.0 

4.6891 

935.4 

1.126 

600.0 

4.6957 

934.8 

1.110 

625.0 

4.7020 

934.2 

1.096 

650.0 

4.7079 

933.5 

1.084 

675.0 

4.7135 

932.9 

1.073 

700.0 

4.7188 

932.3 

1.064 

725.0 

4.7239 

931.6 

1.055 

750.0 

4.7288 

931.0 

1.047 

775.0 

4.7335 

930.4 

1.040 

800.0 

4.7380 

929.8 

1.034 

825.0 

4.7423 

929.1 

1.029 

850.0 

4.7465 

928.5 

1.024 

875.0 

4.7506 

927.9 

1.020 

900.0 

4.7546 

927.2 

1.016 

925.0 

4.7585 

926.6 

1.012 

950.0 

4.7623 

926.0 

1.009 

975.0 

4.7661 

925.4 

1.006 

1000.0 

4.7698 

925.1 

1.003 
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APPENDIX  E;  PROGRAMS 


C 

c 

Q  ***************************+*++}+*++*+******* ************* 

c 

C  THE  COMPRESSIELE  MODEL 
C 

c  *********  *4********** ******** ***************************** 

c 

c 

C  THIS  PROGBAM  CALCULATES  SPELECMICEOCLIMATIC  PARAMETERS 
C  GIVEN  CERTAIN  CHARACT ERISTICS  Of  THE  CAVE  AND  THE 
C  EXTERNAL  ATMOSPHERIC  CONDITIONS.  THE  AIR  IS  ASSUMED 
C  TO  BE  COMPRESSIBLE. 

C 

c 

c 

DOUBLE  PRECISION  M,SVP,MTEST,DABS, ESTMTD, ACTUAL 
DIMENSION  T(101)  ,TO  (101)  ,  TV  (101)  ,D(101)  , DO  (101) 
DIMENSION  P(101)  ,PO(101)  ,QSPEC(101) 

DIMENSION  M  (101)  ,  V  (101)  ,  RHMD  (101) 

DIMENSION  FI  (101)  ,FII(101)  ,FIII(101) 

DIMENSION  X  ( 1 0 1 ) 

REAL  L,LENGTH,K,JI,JII,LAPSE,INPUT 
C 
C 

C  VIRT ( A, B)  CALCULATES  THE  VIRTUAL  TEMPERATURE  AT 
C  TEMPERATURE  A  AND  SPECIFIC  HUMIDITY  B 
C 
C 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


VIRT (A,B)=A* (1.0+0.609*8) 


READ  THE  EXTERNAL  ATMOSPHERIC  CONDITIONS 

TEXTL  =  ATMOSPHERIC  TEMPERATURE  AT  LOWER  ENTRANCE 
PEXTL  =  ATMOSPHERIC  PRESSURE  AT  LOWER  ENTRANCE 
PEXTU  =  ATMOSPHERIC  PRESSURE  AT  UPPER  ENTRANCE 
RH  =  RELATIVE  HUMIDITY  AT  LOWER  ENTRANCE  (DECIMAL) 

N  =  NUMBER  OF  DIVISIONS  AT  WHICH  PARAMETRES  ARE 

CALCULATED  (THE  LARGER  N  IS,  THE  GREATER  THE 
ACCURACY;  MAXIMUM  =  100) 


READ  (4,100)  TEXTL , PEXTL, PEXTU, RH, N 
100  FORMAT  (F6.2,T10,F6. 1,T20,F6 .1 ,T30,F4. 1 ,T40,I4) 
C 

c 

C  READ  THE  CHAR ACTERISTICS  OF  THE  CAVE 
C 

C  LENGTH  =  LENGTH  OF  CAVE 

C  FCOEFF  =  COEFFICIENT  OF  FRICTION  (FANNING) 
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C  RADIOS  =  RADIOS  OF  PASSAGE 
C  TWALL  =  HALL  TEMPERATURE  IN  CAVE 
C  THETA  =  ANGLE  OF  PASSAGE 

C  E  =  FRACTION  OF  PASSAGE  WALL  THAT  IS  MOIST  (0.0  OR 

C  1.0) 

C 

C 


READ (4, 101)  LENGTH, FCOEFF, RADI  OS, T WALL, THET A, E 
101  FORMAT  (F7. 1 ,T 10 , F5. 3,T20,F5 .2,T30 ,F6 . 2 ,T40  ,F4. 1 ,T50 
1  ,F3.1) 


C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 


CONSTANTS 
K 

TFREEZ 
EGAS 
RVAP 
CP 
G 
L 


RATIO  OF  SPECIFIC  HEATS  FOR  DRY  AIR 
FREEZING  POINT  OF  HATER 
GAS  CONSTANT  FOR  DRY  AIR 
GAS  CONSTANT  FOR  HATER  VAPOR 
SPECIFIC  HEAT  CONSTANT  PRESSURE 
ACCELERATION  DOE  TO  GRAVITY 
DRY  ADIABATIC  LAPSE  RATE 


K 

TFREEZ 

RGAS 

RVAP 

CP 

G 

L 


1.4 

273.16 

287.04 

461.50 

1004.0 

9.80 

0.00976 


C 

C 

C  PO(1)  AND  TO  (1)  ARE  THE  INITIAL  STAGNATION  PRESSURE  AND 
C  TEMPERATURE 
C 

C  X(I)  IS  THE  DISTANCE  INTO  THE  CAVE,  WHILE  DELTA  IS  THE 
C  INCREMENT  OE  DISTANCE  BETWEEN  CALCULATED  CROSS-SECTIONS 


C 

c 

PEXIL= 100 . *PEXTL 

PEXTU=  100 .*PEXTU 

PO  (1) =EEXTL 

TO  (1)  =TEXTL+TFREEZ 

XN=N 

NN=N+ 1 

DELI A=LENGTH/ (XN) 
X(1)=0.0 
DO  102  1=2, NN 
X  (I)  =X  (1-1)  +DELTA 
102  CONTINUE 

THALL=TWALL+TFREEZ 
SMOOTH= 1.0 
C 
C 
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C  QWALL 
C 

C  QST 
C 

C  QSPEC(I) 

C  T  (I) 

C  TV  (I) 

C  D(I) 
c  w 
c 
c 

QWALL=SAT  (TWALL,PO  (1) ) 

QST=S AT  (TO  (1)  ,PO(1)  ) 

QSPEC  (  1)  =  (QST*RH)  /  (  (QST*  (RH- 1 . 0) )  ♦  1 . 0) 

C 

C 

C  ASSUME  -INITIAL  VELOCITY  =1.0  AND  CALCULATE  INITIAL 
C  CONDITIONS,  ESPECIALLY  THE  INITIAL  MACH  NUMBER  M(1) 

C 

C 

V  (1)  =  1.0 

103  RE AD  (5,999)  INPUT 
999  FORMAT  (E6. 3) 

IF  (INPUT  „NE.  0.0)  V(1)=INPUT 
T  (1 )  =TO  (1)  +  (V  (1)  **2.0)/  (2 . 0*CP) 

TV  (1)  =VIRT  (T  (1)  , QSPEC  (1)) 

M  (  1)  =V  (1)  *((K*RGAS*TV  (1)  )  **  (-0.5)  ) 

FI  (1 )  =  (M  (1 )  **2.  0)  *  ( 1 . 0  +  K*  (M  (1)  **2.0)  )/ 

1  ( 1 . 0-K*  (M  (1 )  **2 . 0)  ) 

FI  I  ( 1 )  =  (M  (1)  **2.0)  *  (2.0*K*(M  (1)  **2.0) )/ 

1  (1 . 0-K*  (M  ( 1)  **2.0)  ) 

Fill  ( 1)  =  (M  ( 1)  **2.0)  *2. 0/(1 .0-K*  (M  (1)  **2 .0)  ) 
P(1)=PO(1)*((1.0  +(K-1.0)*(M(1)**2.0)/2.0)**((K-1.0) 
1  /K)) 

D(1)=P(1)/(RGAS*TV(1)  ) 

W=D (1) *3. 14159*  (RADIUS**2 . 0) *V  (1) 

C 

C 

C  MONITR  =  MONITORING  DEVICE  FOR  HUMIDITY  (IF  0,  AIR 
C  IS  NOT  YET  SATURATED,  IF  1 ,  IT  IS) 

C 

c 

MONITR  =  0 


=  SPECIFIC  HUMIDITY  AT  SATURATION  FOR  THE 
WALL  TEMPERATURE  AT  ATMOSPHERIC  PRESSURE 
=  SPECIFIC  HUMIDITY  AT  SATURATION  FOR 
INITIAL  ATMOSPHERIC  CONDITIONS 
=  SPECIFIC  HUMIDITY  AT  DISTANCE  X (I)  INTO  CAVE 
=  TEMPERATURE  AT  X  (I) 

=  VIRTUAL  TEMPERATURE  AT  X (I) 

=  AIR  DENSITY  AT  X (I) 

=  MASS  RATE  OF  FLOW 


C 

C 

C 

C 

c 

c 


SMOOTH 

XO 


=  COEFFICIENT  OF  FRICTION  EOR  SMOOTH  PASSAGE 
=  RELAXATION  LENGTH 


IF (FCOEFF.EQ. SMOOTH) FCOEFF=2.0 

SMOOTHS.  25*  (ABS(0. 8685 9*AIOG (  (35886. 1 ADIOS)  / 

1  (1.964*ALOG (35886. 12*H/HADI0S-3. 82 15) ) ) ** 

2  (-2.0)) 

IF  (PCO EPF. GE.  1.000)  FCOEFF  =  SMOOTH 
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X0=(  (SMCOTH/FCOEFF) ♦*0.5) *169.86*  (RADI OS* *0 . 8) * 

1  ( W**0 . 2) 

C 

C 

C  BEGIN  DO  LOOP  TO  CALCULATE  SUBSEQUENT  CROSS-SECTIONS 

C 

C 

DO  110  1=2, NN 
C 
C 

C  IF  THE  TEMPERATURE  IS  BELOW  FREEZING,  H  IS  THE  LATENT 
C  HEAT  OF  SUBLIMATION.  IF  IT  IS  ABOVE  FREEZING,  H  IS 
C  THE  LATENT  HEAT  OF  VAPORISATION 
C 
C 

IF  (T  (1-1)  .LT.TFBEEZ)  H=2.50E6 
IF  (T  (1-1)  .  GE.  TFREEZ)  H=2 . 83E6 
C 
C 

C  THE  TEMPERATURE  IS  NOW  CALCULATED.  THE  PARTICULAR 
C  FORMULA  USED  DEPENDS  ON  WHETHER  OR  NOT  THE  AIR  IS 
C  SATURATED 
C 
C 

IF  (MONITR.NE. 0) GO  TO  105 

QSPEC  (I)  =QWAL1+  (QS PEC  (1-QWALL)  *EXP  (-X  (I)/XO) 

IF  (QSPEC  (I)  .GE.SAT  (T  (1-1)  ,PO  (1)  )  )  GO  TO  104 
T (I) =TWALL-L*XO*SIN (THETA*3. 14159/180.0) 

1  +  (T  (1-TW ALL+L*XO*SIN  (THETA*3 .  1 4 159/1  80 . 0)  )  * 

2  EXP  (-X  (I)/XQ)  ♦  (H*E*  (QSPEC  (1-QWALL)  *X  (I)  * 

3  EXP  (-E*X  (I)/XO)  )/(CP*XO) 

GO  TO  107 

104  XSAT=X  (I—  1) 

TSAT=T  (1-1) 

MONITR= 1 

105  QSPEC  (I)=SAT(T  (1-1)  ,PC(1)) 

T (I) =T  (1-1) 

JI=QW ALL*L/  (CP*RVAP* (TWALL**2 . 0) ) 

JII=L/  (RVAP*  (TWALL**2 . 0) ) 

LAPSE=L*XO*SIN(THETA*3. 14159/180.0) 

1 06  TTEST=TWALL-  (  (1 . 0+JI) *LAPSE/ (JI*JII*LAPSE+ 1.0))  + 

1  (TSAT-TKALL+  (  (1.0+JI)  *LAPSE/  (JI*JII*LAPSE+ 

2  1.0))) *EXP  (  (  (ABS  ( JI*JII*LAPSE+ 1 . 0) * *2. 0) * 

3  (XSAT-X  (I)  -  (T  (I-TSAT)  *XO*JI*JII*  (JI*  JII* 

4  IAPSE+1 .0) ) /  (XO*  (1 . 0+JI) ) ) 

IF  (ABS  (T  (I-TTEST)  .LT.0.01) GO  TO  107 
T  (I) =TTEST 
GO  TO  106 
C 

C  CALCULATE  VIRTUAL  TEMPERATURE  AND  MACH  NUMBER 

C 

C 

1  07  TV  (I)  =VIRT  (T  (I)  , QSPEC  (I) ) 
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M  (I)  =  M  (1-1) 

108  FI  (I)  =  (M  (I)  **2.0)  *  (1.0+K*(M  (I)  **2 . 0)  )/ 

1  (1 .  0-K*  (M  (I)  **2. 0)  ) 

FII  (I)  =  (M  (I)*  *2.0)*  (2.0*K*  (M  (I)  **2 . 0)  )  / 

1  (1.0-K*(M(I)**2.0)) 

Fill  (I)  =  (M  (I)  **2.0)  *2.0/  (1. 0-K*  (M  (I)  **2.0)  ) 

MTEST=  (  (M  (I-1)**2.0)  +  ((  (FI(I)/TV(I)  )  +FI  (I-1)/TV  (1-1)  ) 

1  *0.5*  (TV  (I-TV  (1-1)  )  )  4-  (FCOEFF*  (FII  (I)  +FII  (1-1)  ) 

2  *  (X  (I-X  (1-1)  )/(2. 0*RADIUS)  )  *  (G*SIN  (THETA* 

3  2.0*3.14159/360.0)*  (  (Fill  (I) /TV  (I) )  + 

4  Fill  (I-  1)  /TV  (1-1)  )  *  (X  (I-X  (1-1)  )/(2.  0* 

5  RGAS)  )  )  **  (0 . 5) 

IF  (DABS  (M  (I-MTEST)  .IT.  (1  .E-6)  )  GO  TO  109 
M (I) =MTEST 
GO  TO  108 

109  P  (I)  =P  (1)  *  (M  (1)/M  (I)  )  *  (  (TV  (I) /TV  (1)  )  **0.5) 

1  10  CONTINUE 

WRITE  (6 ,9000)  V(1),P(N+1) 

9000  FORMAT  (1X,F7.4,5X,F9.2) 

C 

c 

C  NOW  ITERATE  FOR  INITIAL  VELOCITY 

C 

C 

ESTMTD=DBLE(P(N+1)  ) 

ACTUAL=EBLE (PEXTU) 

VTEST= V  (1 )  *  (  (ESTMTD/ACTUAL)  **4 . 0) 

IF (ABS  (V  (1-VTEST)  .LT. 0.0001) GO  TO  111 
V  (1)  =VTEST 
GO  TO  103 
C 

c 

C  CALCULATE  ALL  REMAINING  PABAMETRES 

C 

C 

111  RHMD  (1)=QSPEC  (1)  *  (SAT  (T  (1 )  ,  PO  ( 1 )  -  1 . 0)  / 

1  (SAT  (T  (1 )  ,PO(1)  )*  (QSPEC  (1-1.0)) 

DO  112  1=2, NN 

112  RHMD  (I)  =QSPEC  (I)  *  (SAT  (T  (1-1 )  ,PO  (1)  -  1 . 0)  / 

1  (SAT  (T  (1-1)  ,PO(1))*  (QSPEC  <1-1.0)) 

DO  113  1=1, NN 

PO  (I)  =P  (I)  *((1.0+  (K-1 . 0)  *  (M  (I)  **2 . 0) /2 •  0)  ** 

1  (K/  (K-1  .0) )  ) 

D  (I)=P  (I)  /  (RG  AS*T  V  (I)  ) 

DO  (I)  =D  (I)  *  (  ( 1  •  0*  (K-1 . 0)  *  (M  (I)  **2 . 0)  /2 . 0)  **  (1.0/ 

1  (K- 1.0))) 

TO  (I)  =T  (I)  *  (1.0+  (K—  1-0)  *(M  (I)**2.0)/2.0-TFREEZ 
T  (I)  =T  (I-TFREEZ 
TV  (I) =TV (I-TFREEZ 
PO  (I)  =  •  01  *PO  (I) 

P  (I)  =  .01*P(I) 

V(I)=V(1)  *D(1)/D(I) 

113  CONTINUE 

T WALL =TWALL-T FREE Z 
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PEXTU=. 01*PEXTU 
PEXTL= . 01 *PEXTL 
WRITE  (7,114) 

114  FORMAT  ('«*•) 

WRITE  (7,  115)  TEXTL 

115  FORMAT  (IX,1 TEMPERATURE  AT  LOWER  ENTRANCE1  ,T40, 

1  F8.3,T50, •DEG  C. « ) 

WRITE  (7,116)  PEXTL 

116  FORMAT  (IX,1 PRESSURE  AT  LOWER  ENTRANCE' ,T40, 

1  F8. 3 ,T50, 'MB') 

WRITE  (7, 117)  PEXTU 

117  FORMAT  (IX,' PRESSURE  AT  UPPER  ENTRANCE' ,T40, 

1  E 8. 3,T50, ' MB' ) 

WRITE  (7,118)  RH 

118  FORMAT  (IX, 'RELATIVE  HUMIDITY  AT  LOWER  ENTRANCE', 

1  T40,F8 . 3) 

WRITE  (7,119)  N 

119  FORMAT  (IX, 'NUMBER  OF  DIVISIONS  OF  PASSAGE ' ,T40 ,14) 
WRITE  (7,120) 

120  FORMAT  («—*•) 

WRITE  (7,121)  LENGTH 

121  FORMAT  (IX, 'LENGTH  OF  P ASS AGE' , T40 , F8 . 3 , T5 0 , ' M« ) 

WRITE  (7,122)  RADIUS 

122  FORMAT  (IX, 'RADIUS  OF  PASSAGE* ,T40,F8. 3, T50, ' M' ) 

WRITE  (7,123)  THETA 

123  FORMAT  (IX,' ANGLE  OF  P ASSAGE ' , T4 0,F8 . 3 ,T50 , ' DEG ' ) 
WRITE  (7,124)  TWALL 

124  FORMAT  (IX, 'WAIL  TEMPERATURE ' ,T40, F8 . 3, T50 ,' DEG  C.«) 
WRITE  (7,125)  FCOEFF 

125  FORMAT  (IX, 'FRICTION  FACTOR' , 140, F8. 3) 

WRITE  (7,126)  E 

126  FORMAT  ( IX, 'FRACTION  OF  WALL  THAT  IS  MOIST «  ,T40 ,F8 . 3) 
WRITE  (7,127) 

127  FORMAT  ( '-' ) 

WRITE  (7,128) 

128  FORMAT  (1 X ,' DISTANCE ' ,T1 5, 'TEMPERATURE' , T30 , 

1  'PRESSURE' ,T45, 'DENSITY') 

WRITE  (7,129) 

129  FORMAT  (IX, •  (METRES)  • ,T1 5, '  (DEG  CELSIUS)  ', T30, 

1  «  (MILLIBARS)  • ,T45,«  (KG/CUBIC  M) • ) 

WRITE  (7,130) 

130  FORMAT  (.'  —  ') 

DO  132  1=1, NN 

WRITE  (7,131)  X  (I)  ,T  (I)  ,  P  (I)  *  D  (I) 

131  FORMAT  (IX ,F7 • 1 , T15 ,F6. 2 ,T30 ,F6 . 1,T45,F6.4) 

132  CONTINUE 
WRITE  (7,133) 

133  FORMAT  (-'•') 

WRITE  (7 ,134) 

134  FORMAT  (IX ,T1 5 ,' STAGNATION • , T30 , ' STAGNATIO N • ,T45 , 

1  ' STAGNATION ' ) 

WRITE  (7 ,135) 

135  FORMAT  (IX,' DISTANCE' ,T1 5, 'TEMPERATURE' ,T30 , 

1  • PRESSURE' ,T45, ' DENSITY') 


•'  , 

>  . 1  k 

' 

i 

t 

♦ 

. 

• 

•  .. 

* 

.  V 

/ 

* 

■ 

1 

f  1  1 

* 

4 

■ 

P 

> 

•’  *  •  1 

1 

- 

’ 

. 

• 

1 

• 

. 

.  f 

i 

i  ■  '  > 

.  , 

* 

I 

r 

i 

’  ■ 

■ 

! 

♦ 

, 

1 

<  i 

. 

r , 

' 

^  V 

i 

>  V 

« 

«  . 

122 


WRITE  (7,136) 

136  FORMAT  ( IX, •  (METRES)  •  ,T15, •  (DEG  CELSIUS)  • ,T30 , 

1  *  (MILLIEARS)  • ,T45, ■  (KG/CUBIC  M)  '  ) 

WRITE  (7,137) 

137  FORMAT  (  *—  • ) 

DO  139  1=1, NN 

WRITE  (7,138)  X  (I)  ,TO(I)  ,PO  (I)  ,DO(I) 

138  FORMAT  (1X,F7.  1,T15,F6.2 ,T30,F6 . 1,T45,F6. 4) 

139  CONTINUE 
WRITE  (7,140) 

140  FORMAT  ( '-  • ) 

WRITE  (7,141) 

141  FORMAT  ( 1 X, * DI STANCE ',T 15, 'VELOCITY* ,T30,' MACH* , 

1  145,' RELATIVE') 

WRITE  (7,142) 

142  FORMAT  (IX,*  (METRES)  ' ,T15, •  (M/SEC)  • , T30, * N UMBER* , 

1  T45, ' HUMIDITY*) 

WRITE  (7,143) 

143  FORMAT  <•-') 

DO  145  1=1, NN 

WRITE  (7,144)  X  (I)  ,V  (I)  ,M(I)  ,  RHMD  (I) 

144  FORMAT  (IX ,F7 . 1 ,T15 ,F7. 4 ,T30 ,F8 - 6,T45, F5 . 3) 

145  CONTINUE 
WRITE  (7,146) 

146  FORMAT  ( •-  • ) 

STOP 

END 

C 

C 

C  SAT  (A, B)  CALCULATES  THE  SPECIFIC  HUMIDITY  AT  SATURATION 
C  FOR  TEMPERATURE  A  AND  PRESSURE  B 
C 

C  SVP  IS  THE  SATURATION  VAPOR  PRESSURE  AT  TEMPERATURE  T 

C 

C 

FUNCTION  SAT  ( A,  B) 

IF  (A. LT.  273.  16)  GO  TO  1000 

SVP  =10.0**  (-7. 90298*  ((373.  16/A-1 .0) 

1  +  5.02808*ALOG10  (373.  16/A) 

2  -  (1.3816E-7)  *  (10.0**  (11.344*  (1.0-  (A/373.  16) )  ) 

3  -1.0)  +  (8. 1328E-3)*  ((10.0**  (-3.49149* 

4  i  (373.  16/ A- 1 .0)  )-1  .0)  +ALOG10  (1013.  24  6) ) 

GO  TO  1001 

1000  SVP  =10.0**  (-9.09718*  (  (273.  16/A-1 .0) 

1  -3.56654*ALOG1 0(273. 16/A) 

2  +0. 876793* (1 .0- (A/273.  16) ) +ALOG10  (6. 1071) ) 

1001  SAT=1 • 0-  ( 1 . 0/  ( (  (0.621 97*S VP)/ (.01*E-SVP))  +1.0) ) 
RETURN 

END 


. 


•  •• ,  • 

.  i 

*.  V 

J. 

- ' \ 

' 

• 

(  « 1  .  dr  x  ,  / 

, 

. 

"  > 

*  ' 

• 

. 

. 

1  .  i 

» 

• 

* 

* 

■ 

k 

• 

* 

•>  < 

,  '  • 

. 

( 

.  i 

. 

V 

> 

123 


C 

C 

C  ****** ** * ** *  ******* ********** ***  ************ ************** 
C 

C  THE  INCOMPRESSIBLE  MODEL  -  CONSTRAINED 
C 

C  ********************************************************** 

c 

c 

C  THIS  PROGRAM  CALCULATES  SPELECMICBOCLIMATIC  PABAMETEES 
C  GIVEN  CEETAIN  CHABACT EBISTICS  Of  THE  CAVE  AND  THE 
C  EXTERNAL  ATMOSPHERIC  CONDITIONS.  THE  AIR  IS  ASSUMED 
C  TO  BE  INCOMPRESSIBLE. 

C 

C 

C 

DIMENSION  V  (101)  ,£(101)  ,£0(101) 

DIMENSION  T  (101)  ,D  (101)  ,TV  (101) 

DIMENSION  QSPEC  (101)  ,RHMD  (101) 

DIMENSION  X(101) 

REAL  L, LENGTH, JI, JII,INTA, I NTB, LAPSE, INPUT 
C 
C 

C  VIRT  (A , E)  CALCULATES  THE  VIRTUAL  TEMPERATURE  AT 
C  TEMPERATURE  A  AND  SPECIFIC  HUMIDITY  B 


C 

c 

VIRT  ( A  ,E)  =A*(1  .0*0. 609* E) 


C 

C 

C  READ  THE  EXTERNAL  ATMOSPHERIC  CONDITIONS 
C 

C  TEXTL  =  ATMOSPHERIC  TEMPERATURE  AT  LOWER  ENTRANCE 
C  PEXTL  =  ATMOSPHERIC  PRESSURE  AT  LOWER  ENTRANCE 

C  PEXTU  =  ATMOSPHERIC  PRESSURE  AT  UPPER  ENTRANCE 

C  RH  =  RELATIVE  HUMIDITY  AT  LOWER  ENTRANCE 

C  N  =  NUMBER  OF  DIVISIONS  AT  WHICH  PARAMETBES  ARE 

C  CALCULATED  (THE  LARGER  N  IS,  THE  GREATER  THE 

C  ACCURACY;  MAXIMUM  =  100) 

C 

c 

READ  (4,100)  TEXTL, PEXTL, PEXTU, RH, N 
100  FORMAT  (F6 .2,T  10, F6. 1 ,T20,F6 . 1 ,T30,F4 . 1 ,T4  0,I4) 


READ  THE  CHAR ACTERISTICS  OF  THE  CAVE 


C 

c 
c 
c 

C  LENGTH 
C  FCOEFF 
C  RADIUS 
C  TWALL 
C  THETA 
C  E 
C 


=  LENGTH  OF  CAVE 

=  COEFFICIENT  OF  FRICTION  (FANNING) 

=  RADIUS  OF  PASSAGE 
=  WAIL  TEMPERATURE  IN  CAVE 
=  ANGLE  OF  PASSAGE 

=  FRACTION  OF  PASSAGE  WALL  THAT  IS  MOIST  (0.0  OR 

1.0) 


•  . 
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C 

C 


C 

C 

C 

c 


c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


READ  (4/101)  LENGTH , FCOEFF /RADIUS, T WALL, THETA/ E 
101  FORMAT  (F7. 1,T  10,F5.3,T20,F5 .2,T30,F6 .2,T4  0,F4. 1 ,T50 
1  ,F3.1) 


CONSTANTS 


TFREEZ  =  FREEZING  POINT  OF  HATER 

RGAS  =  GAS  CONSTANT  FOR  DRY  AIR 

RVAP  =  GAS  CONSTANT  FOR  HATER  VAPOR 

CP  =  SPECIFIC  HEAT  CONSTANT  PRESSURE 

G  =  ACCELERATION  DUE  TO  GRAVITY 

L  =  DRY  ADIABATIC  LAESE  RATE 


TFREEZ 

RGAS 

RVAP 

CP 

G 

L 


273.  16 

287.04 

461.50 

1004.64 

9.80 

0.00976 


PO(1)  IS  THE  INITIAL  STAGNATION  PRESSURE 

\ 

X(I)  IS  THE  DISTANCE  INTO  THE  CAVE,  WHILE  DELTA  IS  THE 
INCREMENT  OF  DISTANCE  BETWEEN  CALCULATED  CROSS-SECTIONS 


PEXTL=  100.*PEXTL 
PEXTU=100.*PEXTU 
PO  (1)  =EEXTL 
X  N=  N 
NN  =  N  + 1 

DELTA=LENGTH/(XN) 

X (1)  =0.0 
DO  102  1=2, NN 
X  (I)=X  <I-1)  +DELTA 
102  CONTINUE 

THAIL=T  WALL+TFREEZ 
SMCOTH= 1.0 


QWALL 

QSPEC  (I) 
QST 

T(I) 

TV  (I) 

D  (I) 

W 


SPECIFIC  HUMIDITY  AT  SATURATION  FOR  THE 
WALL  TEMPERATURE 

SPECIFIC  HUMIDITY  AT  DISTANCE  X (I)  INTO  CAVE 
SPECIFIC  HUMIDITY  AT  SATURATION  FOR 
INITIAL  ATMOSPHERIC  CONDITIONS 
TEMPERATURE  AT  X  (I) 

VIRTUAL  TEMPERATURE  AT  X (I) 

AIR  DENSITY  AT  X  (I) 

MASS  RATE  OF  FLOW 
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C 

T  (1 )  =TEXTL+TFREEZ 
QWALL=SAT  (TWALL,PO  (1)  ) 

QST-SAT  (T  (1)  ,  PO  (1 )  ) 

QSPEC  ( 1)  =  (QST*RH)  /  (  (QST*  (RH-1  . 0)  )  +1 .0) 

c 

C 

C  CALCULATE  THE  VIRTUAL  TEMPERATURE  AND  THE  AIR 
C  DENSITY  AT  THE  LOWER  ENTRANCE 
C 
C 

TV  (1 )  =  VIRT  (T  ( 1)  , QSPEC  ( 1 )  ) 

D(1)=PG  (1)/(RGAS*TV  (1)) 

C 

C 

C  ASSUME  INITIAL  VELOCITY  =  1.0  AND  CALCULATE  INITIAL 
C  PRESSURE 
C 

c 

V(1)  =  1.0 

103  READ  (5  *999)  INPUT 
999  FORMAT  (F6. 3) 

IF  (INPUT  .NE.  0.0) V  (1) =INPUT 
P(1)  =PO  (1-0. 5*D  (1)  *  (V  (1 )  **2 . 0) 

W=D  (1) *3. 14159* (RADIUS* *2 . 0) *V(1) 

C  ' 

C 

C  MONITR  =  MONITORING  DEVICE  FCR  HUMIDITY  (IF  0,  AIR 
C  IS  NOT  YET  SATURATED ,  IF  1 ,  IT  IS) 

C 

C 

MONITR  =  0 
C 

c 

C  SMOOTH  =  COEFFICIENT  OF  FRICTION  FOR  SMOOTH  PASSAGE 

C  XO  =  RELAXATION  LENGTH 

C 
C 

IF (ECOEFF.EQ. SMOOTH) FCOEFF=2.0 

SMGOTH=0. 25*  (ABS  (0 . 86859*ALOG (  (35886 . 1 2* W/fi ADIUS)  / 

1  (1 .964*ALOG  (35886 .  12*W/RADIUS-3. 82 1 5) ) ) ** 

2  (-2.0)  ) 

IF  (ECOEFF.GE.  1.000) ECOEEF=SMCOTH 

XO=  (  (SMCOTH/FCOEFF)  **0.5)  *169.86*  (RADIUS**0. 8)  * 

1  (W**0. 2) 

C 

C 

C  BEGIN  DC  LOOP  TO  CALCULATE  SUBSEQUENT  CROSS-S ECTIONS 

C 

C 

DO  110  1=2, NN 


IF  THE  TEMPERATURE  IS  BELCW  FREEZING,  H  IS  THE  LATENT 


•  ,  «  i4*2i  =  U  t  > 
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C  HEAT  OF  SUELIMATION.  IF  IT  IS  ABOVE  FREEZING ,  H  IS 
C  THE  LATENT  HEAT  OF  VAPORISATION 
C 
C 

IF (T  (1-1) -LT.TFREEZ) H=2.50E6 
IF  (T  (1-1)  .GE.TFREEZ)  H=2.83E6 
C 
C 

C  THE  TEMPERATURE  IS  NOW  CALCULATED.  THE  PARTICULAR 
C  FORMULA  USED  DEPENDS  ON  WHETHER  OR  NOT  THE  AIR  IS 
C  SATURATED 
C 
C 

IF  (MONITR.NE. 0)  GO  TO  105 

QSPEC  (I)  =QWALL+  (QSPEC  (1-QWALL)  *EXP  (-X  (I)  /XO) 

IF  (QSPEC  (I)  .GE.SAT  (T  (1-1)  ,PO  (1)  )  )  GOTO  104 
T  (I) =TWALL-L*XO*SIN  (THETA*3 . 14159/180.0) 

1  +  (T  ( 1-T WALL+  L*XO*SIN  (THETA*3. 14159/1 80.0) ) * 

2  EXP (— X (I)/XO)  +  (H*E* (CS PEC (1-QWALL)  *X (I) * 

3  EXP  (-E*X  (I)  /XO)  )  /  (CP*XO) 

GO  TO  107 

104  XSAT=X  (I-  1) 

TSAT=T  (1-1) 

MCNITR=1 

105  T  (I)  =T  (1-1) 

QSPEC  (I)  =SAT  (T  (I-  1)  ,PO  (1)  ) 
JI=QWAIL*L/(CP*RVAP*(TWALL**2.0)  ) 

JII=L/  (BVAP* (TW ALI**2 . 0)  ) 

LAPSE=L*XO*SIN(THETA*3. 14159/180.0) 

1 06  TTEST=T WALL- ( (1 . 0+JI) *LAPSE/ (JI*JII*LAPSE+1 . 0)  )  + 

1  (TSAT-TWALL  +  (  (1 .0+JI) ELAPSE/  (JI*JII*LAPSE+ 

2  1.0))) *EXP  (  (  (ABS ( JI* JII*LAPSE+ 1 . 0) **2. 0)  * 

3  (XSAT-X(I)-  (T  (I-TSAT) *XO*JI*JII* (JI*JII* 

4  LAPSE+1 .0) ) /  (XO*  ( 1 - 0+JI) ) ) 

IF  (AES  (T  (I-TTEST)  .LT.0.C1)GC  TO  107 
T  (I)  =1TEST 
GO  TO  106 
C 
C 

C  CALCULATE  VIRTUAL  TEMPERATURE 

C 

C 

107  TV  (I)  =VIRT  (T  (I)  #QSPEC  (I)  ) 


CALCULATE  DENSITY 


1 

2 

3 

1 


D  (I)=D  (1)  *  (TV  (1)  /TV  (I) )  *EXP  (  (-G*SIN  (THETA 
*3.  14 159/1  80.0)  /EGAS)  *(  (1.0/TV  (I)  )  + 

(1 .0/TV  (1-1)  )  )  *  (  (X  (I-X  (1-1)  )/2. 0)  ) 

*(P(I-1)/PO))  „  _ 

P  (I)  =P  (I-  1- D  (  1)  *  ( V  ( 1)  **2)  *  (  (D  (1 )  /D  (I)  - 1 . 0- 

G*SIN  (THETA*3. 1  41  59/1  80  •  0)  *  (D  (I)  *D  (I“1) ) 
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2  0.5*  (X  (I-X  (I-I)-FCOEFF*  (D  (1)  **2)*(V(1) 

3  ** 2)  *  (  (  1 . 0/D  (I)  )  +  (  1 . 0/D  (I-  1 )  )  )  *  ( X  (I-  X  (I- 1 )  )  / 

4  (2 . 0*R ADIUS) 

110  CONTINUE 
C 

C 

C  NOW  ITERATE  EOR  INITIAL  VELOCITY 

C 

C 

INTA=D  (1)  +  (4.  0*D  (2)  )  +D(N  +  1) 

INTB=  ( 1  . 0/D  ( 1 ) )  +  (4 . 0/D  (2)  )  +  (1.0/D  (N  +  1)  ) 

DO  3000  1=4, N ,  2 

INTB=INTB+  (2. 0/D  (1-1)  )  +  (4.0/D  (I)  ) 

INTA=INTA*  (2.0*D  (1-1)  )  +  (4. 0*D  (I) ) 

3000  CONTINUE 

P(N+1)=PEXTL-D(1)*((D(1)/D(N+1)-0.5)*(V(1)** 

1  2-G*SIN  (THETA*3.  14159/180.0) *INTA* 

2  DELTA/3. 0-FCOEFF*(D  ( 1)  **2)  *  (V  (1)  **2)  * 

3  INTB*DELTA/(3.0*RADIUS) 

VTEST=  V  (1)*  (  (P  (N+1)/PEXTU)**4) 

WRITE  (6,2000)  V(1)  ,P(N+1) 

2000  FORMAT  (1X,F7. 4, 5X ,F8. 2) 

IF (ABS  ( V  ( 1— VTEST)  . LT. 0. 0001 ) GO  TO  111 
V (1)  = VTEST 
GO  TO  103 
C 
C 

C  CALCULATE  ALL  REMAINING  PARAMETRES 

C 

C 

111  RHMD  (1 )  =QSPEC  (1)*  (SAT(T(1)  ,PO(1)-1.0)/ 

1  (SAT  (T  (1)  ,PO(1)  )  *  (QSPEC  (1-1.0)  ) 

DO  112  1=2, NN 

1  12  RHMD  (I)  =QSPEC  (I)  *  (SAT  (T  (I- 1 )  ,  PO  (1 )  -  1 . 0)  / 

1  (SAT  (T  (I- 1 )  ,PO ( 1 ) ) *  (QSPEC  (1-1.0)) 

DO  113  1=1, NN 
V(I)=V(1)*D(1)/D  (I) 

PO  (I)  =P  (I)  +0. 5*D  (I)  *  (V  (I)  **2.0) 

T  (I)  =T  (I-TFREEZ 
TV  (I)  =T V  (I-TFREEZ 
PO(I)=.01*PO(I) 

P  (I)  =  .  0 1*P  (I) 

113  CONTINUE 

T W ALL =T WALL-T FREE Z 
PEXTU= .01 *PEXTU 
PEXTL=. 01*PEXTL 
WRITE  (7,114) 

114  FORMAT  ( •-  •) 

WRITE  (7, 115)  TEXTL 

115  FORMAT  (IX, 'TEMPERATURE  AT  LOWER  ENTRANCE' ,T40, 

1  F8. 3,T50, ' DEG  C.') 

WRITE  (7,116)  PEXTL 

116  FORMAT  (IX, 'PRESSURE  AT  LOWER  ENTRANCE' , T4 0 , 

1  F8.3,T50, 'MB') 
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WRITE(7,117)  PEXTU 

117  FORMAT  (IX, 'PRESSURE  AT  UPPER  ENTRANCE' ,T40, 

1  F8.3,T50,'MB') 

WRITE  (7/118)  RH 

118  FORMAT  (IX, 'RELATIVE  HUMIDITY  AT  LOWER  ENTRANCE' , 

1  T40,F8.3) 

WRITE  (7,119)  N 

119  FORMAT  (IX, 'NUMBER  OF  DIVISIONS  OF  PASSAGE ' ,T40 ,14) 
WRITE  (7,120) 

120  FORMAT  ('-•) 

WRITE  (7,121)  LENGTH 

121  FORMAT  (IX, 'LENGTH  OF  PASSAGE' ,T40, F8. 3, T50, ' M' ) 

WRITE  (7,122)  RADIUS 

122  FORMAT  (IX, 'RADIUS  OF  PASSAGE' ,T40,F8. 3, T50 ,' M' ) 

WRITE  (7,123)  THETA 

123  FORMAT  (IX,' ANGLE  OF  PASSAGE ' ,T4 0,F8 .3 ,T50 ,' DEG' ) 
WRITE  (7,124)  TWALL 

124  FORMAT  (IX,' WAIL  TEMPERATURE ' ,T40,F8 . 3, T50 ,' DEG  C.') 
WRITE  (7, 125)  FCOEFF 

125  FORMAT  (IX, 'FRICTION  FACTOR' ,T40,F8. 3) 

WRITE  (7,126)  E 

126  FORMAT  (IX, 'FRACTION  OF  WALL  THAT  IS  MOIST • ,T40 ,F8. 3) 
WRITE  (7,127) 

127  FORMAT (•-•) 

WRITE  (7,128) 

128  FORMAT  (IX, • DISTANCE' ,T15, 'TEMPERATURE' , T30, 

1  'PRESSURE' ,T45, 'DENSITY') 

WRITE  (7,129) 

129  FORMAT  (IX,'  (METRES)  « ,T15, •  (DEG  CELSIUS)  ', T30, 

1  •  (MILLIBARS)  » ,T45,«  (KG/CUBIC  M) • ) 

WRITE  (7,130) 

130  FORMAT  ( '- ' ) 

DO  132  1=1, NN 

WRITE  (7,131)  X  (I)  ,T  (I)  ,  P  (I)  ,D  (I) 

131  FORMAT  (1X,F7. 1,T15,F6.2,T30,F6. 1,T45,F6.4) 

132  CONTINUE 
WRITE  (7,133) 

133  FORMAT  (•-•) 

WRITE  (7,134) 

134  FORMAT  (IX, T30, 'STAGNATION') 

WRITE  (7,135) 

135  FORMAT  (1 X ,' DISTANCE ' ,T1 5, ' VELOCITY • ,T30, 

1  'PRESSURE' ,T45, 'RELATIVE') 

WRITE  (7,136) 

136  FORMAT  (IX,'  (METRES) • ,T15,  '  (M/SEC)  ' ,T30, 

1  «  (MILLIBARS)  '  ,T45,« HUMIDITY') 

WRITE  (7,137) 

137  FORMAT  ('-—') 

DO  139  1=1, NN 

WRITE  (7,138)  X(I)  ,V  (I)  ,PO(I)  ,  RHMD  (I) 

138  FORMAT  (1X,F7. 1 ,T1 5 , F7 . 4,T30 ,F6 . 1 , T45,  F5 . 3) 

139  CONTINUE 
WRITE  (7,140) 

140  FORMAT  ( '- ' ) 
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STOP 

END 

C 

c 

C  S AT ( A, E)  CALCULATES  THE  SPECIFIC  HUMIDITY  AT  SATURATION 
C  FOR  TEMPERATURE  A  AND  PRESSURE  B 
C 

C  S VP  IS  THE  SATURATION  VAPOR  PRESSURE  AT  TEMPERATURE  T 

C 

C 

FUNCTION  SAT  (A,  B) 

IF  (A. LT. 273.  16)  GO  TO  1000 

SVP  = 10. 0** (-7.90298* ( (373. 16/A-1 .0) 

1  +5.02808*ALOG10(373. 16/A) 

2  -  (1.381 6E-7)  *  (10.0**  (11.344*  (1.0-  (A/373. 16) )  ) 

3  -1.0)  +  (8. 1328E-3)  *  ( (10.0**  (-3.4914  9* 

4  (  (373. 16/A-1 .0) ) -1 .0) +ALOG10  (1013.  246) ) 

GO  TO  1001 

1000  SVP  =10.0**  (-9. 0971 8* ((273. 16/A-1 . 0) 

1  — 3 . 5 6 654* A LOG  10(273.  16/A) 

2  +0.876793*  (1 . 0- (A/273 . 16)  )  +ALOG10(6.  1071)  ) 

1001  SAT=1 .0-  (1.0/((  (0.62197*SVP)/(.01*B-SVP) ) +1.0) ) 

RETURN 

END 

C 

C 

c  *******  ****** ************************ ******** ************* 

C 

C  THE  INCOMPRESSIBLE  MODEL  -  UNCONSTRAINED 
C 

C  ***********4********************************************** 

C 

C 

C  THIS  PROGRAM  CALCULATES  SPELEOMICROCLIMATIC  PARAMETRES 
C  GIVEN  CERTAIN  CHARACTERISTICS  OF  THE  CAVE  AND  THE 
C  EXTERNAL  ATMOSPHERIC  CONDITIONS.  THE  AIR  IS  ASSUMED 
C  TO  BE  INCOMPRESSIBLE. 

C 

C 

C 

DIMENSION  V  (101)  ,P  (101)  ,PC  (101) 

DIMENSION  T  (101)  ,  D  (101)  ,TV  (101) 

DIMENSION  QSPEC (101)  ,QSAT  (1  0 1)  ,RHMD  (1 0 1 ) 

DIMENSION  X(101) 

REAL  L, LENGTH, INTA,INTB, INPUT 
C 
C 

C  VIRT  (A , E)  CALCULATES  THE  VIRTUAL  TEMPERATURE  AT 
C  TEMPERATURE  A  AND  SPECIFIC  HUMIDITY  B 
C 

c 

VIRT  (A, B)  =  A*(1.0+0.609*E) 

C 

C 


* 


. 

. . 

'  ■ 

- . ,  . 

> 

. 

.  \  ( 

*8-*l  -  X 

.  •  .. 


(  . 

. 


.  • 

, 


. 


. 


. 

. 

♦  ■  i 

( i  A  ) 


'  ■  >  •.! 


.■»  • 


j 


130 


C  READ  THE 
C 

C  TEXTL  = 
C  PEXTL  = 
C  PEXTU  = 
C  RH 
C  N 
C 

c 

c 

c 


EXTERNAL  ATMOSPHERIC  CONDITIONS 

ATMOSPHERIC  TEMPERATURE  AT  LOWER  ENTRANCE 
ATMOSPHERIC  PRESSURE  AT  LOWER  ENTRANCE 
ATMOSPHERIC  PRESSURE  AT  UPPER  ENTRANCE 
RELATIVE  HUMIDITY  AT  LOWER  ENTRANCE 
NUMBER  OF  DIVISIONS  AT  WHICH  PARAMETRES  ARE 
CALCULATED  (THE  LARGER  N  IS,  THE  GREATER  THE 
ACCURACY;  MAXIMUM  =  100) 


READ (4, 100)  TEXTL, PEXTL, PEXTU, RH,N 
100  FORMAT  (F6.2,T10,F6. 1,T20,F6.1 ,T30,F4. 1,T40,I4) 


C 

C 

C  READ  THE  CHARACTERISTICS  OF  THE  CAVE 
C 

C  LENGTH  =  LENGTH  OF  CAVE 

C  FCOEFF  =  COEFFICIENT  OF  FRICTION  (FANNING) 

C  RADIUS  =  RADIUS  OF  PASSAGE 
C  TWALL  =  WALL  TEMPERATURE  IN  CAVE 
C  THETA  =  ANGLE  OF  PASSAGE 

C  E  =  FRACTION  OF  PASSAGE  WALL  THAT  IS  MOIST  (0.0  OR 

C  1.0) 

c 
c 

READ  (4,101)  LENGTH, FCOEFF, RADI US, T WALL, THET A, E 
101  FORMAT  (F7.1 ,T 10 , F5. 3 ,T20, F5 .2 ,T30 , F6 . 2, T40 , F4. 1 ,T50 
1  ,13.1) 

C 

C 

C  CONSTANTS 
C 

C  TFREEZ  =  FREEZING  POINT  OF  WATER 
C  RGAS  =  GAS  CONSTANT  FOR  DRY  AIR 
C  CP  =  SPECIFIC  HEAT  CONSTANT  PRESSURE 

C  G  =  ACCELERATION  DUE  TO  GRAVITY 

CL  =  DRY  ADIAEATIC  LAPSE  RATE 


C 


C 


TFREEZ  = 

273.  16 

RGAS 

287. 04 

RVAP 

461.50 

CP 

1004.64 

G 

9.80 

L 

0.00976 

c 

c 

c 

PO(1)  IS  THE 

INITIAL 

c 

c 

X  (I)  IS  THE 

DISTANCE 

c 

INCREMENT  OF 

DISTANCE 

C 

c 


STAGNATION  PRESSURE 

INTO  THE  CAVE,  WHILE  DELTA  IS  THE 
BETWEEN  CALCULATED  CROSS-SECTIONS 
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PEXTL=100.*PEXTL 
PEXTU=100.*PEXTU 
PO  (1)  =PEXTL 
XN=N 
NN  =  N  + 1 

DELTA=LENGTH/  (XN) 

X(1)=0.0 
DO  102  1=2 , NN 
X  (I)=X  (1-1)  +DELTA 

102  CONTINUE 
T WALL =T HALL +T FREE Z 
SMCCTH= 1.0 

C 

c 

C  QWALL 
C 

C  QSPEC(I) 

C  QST 
C 

C  T  (I) 

C  TV  (I) 

C  D  (I) 

C  W 

c 

T  (1)  =TEXTL+TFREEZ 
QWALL=SAT  (TWALL,PO  (1)  ) 

QST=S  AT  (T  ( 1)  ,PO(1)) 

QSPEC  (1)=(QST*RH)/((QST*<RH-1.0))  +1.0) 

C 

C 

C  CALCULATE  THE  VIRTUAL  TEMPERATURE  AMD  THE  AIR 
C  DENSITY  AT  THE  LOWER  ENTRANCE 
C 
c 

TV  (1)  =VIRT  (T  (  1)  ,  QSPEC  (1)  ) 

D  (1)  =PC  (1)/  (RGAS*TV  (1)  ) 

C 

C 

C  ASSUME  INITIAL  VELOCITY  =1.0  AND  CALCULATE  INITIAL 
C  PRESSURE 
C 
C 

V  (1)  =1 .0 

103  READ  (5,999)  INPUT 
999  FORMAT  (E6. 3) 

IF  (INPUT  .NE.  0.0) V (1) =INPUT 
P  (1)=FC  (1-0.5*D(1)  *(V(1)**2.0) 

QS  AT  (1)  =SAT  (T  (1)  ,PO  (1)  ) 

W=D  (1) *3. 14159* (RADI US* *2.0)  *V(1) 

C 

C 

C  SMOOTH  =  COEFFICIENT  CF  FRICTION  FOR  SMOOTH  PASSAGE 
C  XO  =  RELAXATION  LENGTH 

C 


=  SPECIFIC  HUMIDITY  AT  SATURATION  FOR  THE 
WALL  TEMPERATURE  AT  ATMOSPHERIC  PRESSURE 
=  SPECIFIC  HUMIDITY  AT  DISTANCE  X (I)  INTO  CAVE 
=  SPECIFIC  HUMIDITY  AT  SATURATION  FOR 
INITIAL  ATMOSPHERIC  CONDITIONS 
=  TEMPERATURE  AT  X  (I) 

=  VIRTUAL  TEMPERATURE  AT  X  (I) 

=  AIR  DENSITY  AT  X (I) 

=  MASS  RATE  OF  FLOW 
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C 

IF (ECOEFF.EQ. SMOOTH) FCOEFF=2.0 

SMOOT H=0. 25* (ABS  (0 . 86 859* ALOG < (35886 . 1 2*W/2ADIUS)  / 

1  (1.964*ALOG  (35886 . 12*W/RADIUS-3. 8215) ) ) ** 

2  (-2.0)) 

IF  (ECO EEF.GE.  1.00  0) FCOEFF=S MCOTH 

XO=  (  (SMCOTH/FCOEFF) **0.5) *169.86* (RADIUS**0. 8) * 

1  (W**0 . 2) 

C 

c 

C  BEGIN  DO  LOOP  TO  CALCULATE  SUBSEQUENT  CROSS-SECTIONS 

C 

C 

DO  110  1=2, NN 
C 
C 

C  IF  THE  TEMPERATURE  IS  BELCH  FREEZING,  H  IS  THE  LATENT 
C  HEAT  OF  SUELIMATION .  IF  IT  IS  ABOVE  FREEZING,  H  IS 
C  THE  LATENT  HEAT  OF  VAPORISATION 
C 
C 

IF  (T  (1-1)  . LT. TFREEZ) H=2.50E6 
IF (T  (I- 1)  . GE. TFREEZ) H=2.83E6 
C 
C 

C  THE  TEMPERATURE  IS  NOW  CALCULATED 
C 

c 

QSPEC  (I)=QWALL+  (QSPEC  (1-QWALL)  *EXP  (-X  (I)  /XO) 

L=0. CCS76 

T  (I) =TWALL-L*XO*SIN  (THETA*3 . 14159/180.0) 

1  + (T ( 1— T WALL  +  L*XO*SIN  (THETA*3 . 14 159/1 80 . 0) ) * 

2  EXP  (-X(I)/XO)  +  (H*E*  (QSPEC  (1-QWALL)  *X  (I)  * 

3  EXP  (-E*X  (I)/XO)  )/(CP*XO) 

C 

C 

C  CALCULATE  VIRTUAL  TEMPERATURE 

C 

C 

107  TV  (I)  =VIRT  (T  (I)  ,QSPEC  (I)  ) 

C 

C 

C  CALCULATE  DENSITY 

C 

C 

D  (I)  =D  (1)  *  (TV  (1 )  /TV  (I)  )  *EXP  (  (-G*SIN  (THETA 

1  *3.141 59/1 80 . 0) /RG AS)  *  (  (1.0/TV  (I)  )  ♦ 

2  ( 1 . 0/T  V  ( I-  1 )  )  )  *  (  (  X  (I- X  (I-  1 ) )  /2 . 0 )  ) 

3  *(P(I-1)/P(1))  „ 

P  (I)  =P  (I-1-D  (1)  *  (V  (1)  **2)  *  (  (D  (1)/D  (I)  -1.0- 

1  G*SIN  (THETA* 3.  14159/180.0)  *  (D  (I)  *D  (I- 1) )  * 

2  0.5*  (X  (I-X  (1-1)  -FCCEFF*  (D  (1)  **2)  *  (V  (  1) 

3  **2)  *  (  (1.0/D  (I) )  ♦  (1 . 0/D  (1-1)  )  )  *  (X  (I-X  (1-1)  )  / 

4  (2.0*RADIUS) 
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110  CONTINUE 
C 

C 

C  NOW  ITEBATE  FOB  INITIAL  VELOCITY 

C 

C 

INTA=D  (1)  +  (4.  0*D  (  2)  ) +D  (N+1) 

INTB=  (1.0/D (1))  +  (4. 0/D  (2)  )  +  (1 . 0/D  (N  +  1)  ) 

DO  3000  1=4, N, 2 

INTB=INTB+  (2.  0/D(I-1)  )  +  (4.0/D  (I)  ) 

INTA=INTA+  (2.  0*D(I-1)  )  +  (4 .0*D  (I)  ) 

3000  CONTINUE 

P  (N  +  1)  =PEXTL-D  (1)  *  (  (D  (1)/D  (N+1)-0.5)*(V(1)  ** 

1  2-G*SIN(THETA*3. 14159/180.0) *INTA* 

2  DELTA/3. 0-FCOEFF*  (D(1)  **2)  *  (V  (1)  **2)  * 

3  INTB*DELTA/ (3.0*  BADI US) 

VTEST- V  (1)  *((F(N+1)/PEXTU)  **4) 

WRITE  (6,2000) V ( 1 )  ,P  (N+1) 

2000  FOEMAT  (1X,F7. 4,5X,F8.2) 

IF  (AES  (V  (1-VIEST)  .LT.O. 0001) GO  TO  111 
V (1) =VTEST 
GO  TO  103 
C 
C 

C  CALCULATE  AIL  REMAINING  PAEAMETEES 

C 

C 

111  RHMD  ( 1 )  =QSPEC  (1 )  *  (SAT(T(1)  ,P(1)-1.0)/ 

1  (SAT  (T  (  1)  ,  P  (1 )  )  *(CSFEC  (1-1.0)) 

DO  112  1=2, NN 

112  RHMD  (I)  =QSPEC  (I)  *  (SAT  (T  (1-1)  ,P  (I- 1)  - 1 . 0)  / 

1  (SAT (T  (1-1)  ,P (1-1) ) *  (QSPEC  (1-1.0) ) 

DO  113  1=1, NN 
V  (I)=V  (1)  *D(1)/D(I) 

PO(I)  =P  (I)  +  0.5*D(I)  *  (V  (I)  **2.0) 

T  (I)  =1  (I-TFREEZ 
TV (I) =TV  (I-TFREEZ 
PO(I)  =.01*PO(I) 

P  (I)  =  .  0 1*P  (I) 

1 13  CONTINUE 

T WALL=T WALL-TFEEEZ 
PEXTU= . 0 1*PEXTU 
PEXTL= .01 *PEXTL 
WRITE  (7,114) 

1  14  FORMAT  (•-') 

WRITE  (7,115)  TEXT L 

115  FORMAT  (IX, 'TEMPERATURE  AT  LOWER  ENTRANCE' ,T40, 

1  F8.3,T50,'DEG  C.') 

WRITE  (7,116)  PEXTL 

116  FORMAT ( IX, 'PRESSURE  AT  LOWER  ENTRANCE' ,T4 0 , 

1  F8.3,T50, 'MB') 

WRITE  (7,117)  PEXTU 

117  FORMAT  (IX, 'PRESSURE  AT  UPPER  ENTRANCE' ,T40, 

1  F8.3,T50, ' ME') 
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WRITE  (7,118)  RH 

118  FORMAT  ( IX ,' RELATIVE  HUMIDITY  AT  LOWER  ENTRANCE* , 

1  140, F8. 3) 

WRITE  (7,119)  N 

119  FORMAT  (IX, "NUMBER  OF  DIVISIGNS  OF  PASSAGE ' ,T40,I4) 
WRITE  (7,120) 

120  FORMAT  ('-•) 

WRITE  (7,121)  LENGTH 

121  FORMAT  (IX, ‘LENGTH  OF  PASS AGE' ,T40 , F8 . 3, T5 0 , * M • ) 

WRITE  (7,122)  RADIUS 

122  FORMAT  (IX, 'RADIUS  OF  PASS AG E« , T40 , F8 . 3 , T5 0 , ' M ' ) 

WRITE  (7,123)  THETA 

123  FORMAT  (IX, 'ANGLE  CF  P ASSAGE ' ,T4 0 , F8 . 3 ,T50 , ' DEG  • ) 
WRITE  (7,124)  TWALL 

124  FORMAT  (IX,' WALL  TEMPERATURE' ,T4 0,F8 .3 ,T50 ,' DEG  C.  •) 
WRITE  (7,125)  FCOEFF 

125  FORMAT  (IX, 'FRICTION  FACTOR' ,T40,F8. 3) 

WRITE  (7,126)  E 

126  FORMAT  (IX, 'FRACTION  OF  WALL  THAT  IS  MOIST •  ,T40 ,F8 . 3) 
WRITE  (7,127) 

127  FORMAT  (•—  '  ) 

WRITE  (7,128) 

128  FORMAT  (IX, • DISTANCE' ,T15, 'TEMPERATURE'  ,T30, 

1  'PRESSURE' ,145, 'DENSITY') 

WRITE  (7,129) 

129  FORMAT  (IX,'  (METRES)  • ,T15, '  (DEG  CELSIUS)  ', T30, 

1  •  (MILLIEARS)  ' ,T45,'  (KG/CUBIC  M)  ' ) 

WRITE  (7,130) 

130  FORMAT  ('-') 

DO  132  1=1, NN 

WRITE  (7,131)  X(I)  ,T  (I)  ,P(I)  ,D(I) 

131  FORMAT  (1X,F7. 1,T15,F6.2  ,T30 ,F6 . 1 ,T45, F6. 4) 

132  CONTINUE 
WRITE  (7,133) 

133  FORMAT  (•*—•) 

WRITE  (7,134) 

134  FORMAT  (IX, T30, 'STAGNATION' ) 

WRITE  (7,135) 

135  FORMAT  ( 1 X ,' DISTANCE ' ,T1 5, 'VELOCITY' ,T30, 

1  'PRESSURE' ,T45, 'RELATIVE' ) 

WRITE  (7,136) 

136  FORMAT  (IX,'  (METRES)  •  ,T15,  •  (M/SEC)  '  ,T30, 

1  •  (MILLIBARS)  ' ,T45 ,' HUMIDITY' ) 

WRITE  (7,137) 

137  FORMAT  (•-') 

DO  139  1=1, NN 

WRITE  (7,138)  X  (I)  ,  V  (I)  ,PO  (I)  ,  RHMD  (I) 

1 38  FORMAT  (1X,F7. 1,T15,F7.4,T30,F6. 1,T45,F5.3) 

139  CONTINUE 
WRITE  (7,140) 

140  FORMAT (•-•) 

STOP 

END 
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